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Report  Title 

Globally  convergent  numerical  methods  for  coefficient  inverse  problems 

ABSTRACT 

Coefficient  Inverse  Problems  (CIPs)  for  Partial  Differential  Equations  (PDEs)  represent  a  very  important  tool  for  such  needs  of  the  Army  as 
imaging  of  unknown  targets  hidden  in  cluttered  heterogeneous  backgrounds.  The  goal  of  this  project  is  the  development  of  globally 
convergent  numerical  methods  for  a  wide  class  of  CIPs.  These  methods  are  tested  on  mathematical  models  of  the  interest  to  the  Army  such 
as  imaging  of  antipersonnel  land  mines  and  targets  on  battlefields  covered  by  smogs  and  flames.  In  our  definition  "global  convergence" 
entails:  (1)  a  rigorous  convergence  analysis  that  does  not  depend  on  the  quality  of  the  initial  guess,  and  (2)  numerical  simulations  that 
confirm  the  advertised  convergence  property. 

A  conventional  way  to  solve  a  CIP  is  via  the  minimization  of  a  least  squares  objective  functional.  This  functional  characterizes  misfit 
between  the  data  and  the  solution  of  that  PDE  with  a  "guess"  for  the  unknown  coefficient.  However,  it  is  well  known  to  researchers 
working  on  computations  of  inverse  problems  that  the  phenomenon  of  multiple  local  minima  of  these  functionals  represents  the  major 
obstacle  for  the  development  of  reliable  numerical  methods  for  multidimensional  CIPs.  This  phenomenon  in  turn  is  caused  by  the  above 
mentioned  non-linearity  and  ill-posedness.  Because  of  local  minima,  one  should  somehow  guess  in  advance  about  a  good  approximation  for 
the  solution.  Without  the  availability  of  a  first  good  guess,  however,  there  is  no  guarantee  that  the  calculated  coefficient  is  indeed  close  to 
the  correct  one.  In  our  tenninology  these  are  locally  convergent  numerical  methods.  In  other  words,  their  convergence  to  the  correct  solution 
can  be  guaranteed  only  if  the  starting  point  is  located  in  a  small  neighborhood  of  this  solution.  Because  of  local  minima,  conventional 
numerical  methods  for  multidimensional  CIPs  are  locally  convergent  ones.  However,  in  many  important  applications  the  first  good  guess  is 
unavailable.  In  particular,  locally  convergent  algorithms  are  \QTR{bf}  {fundamentally  unsatisfactory}  for  the  needs  of  the  Army,  because  an 
accurate  a  priori  knowledge  of  the  properties  of  a  medium  is  rarely  available  in  military  applications.  This  is  because  military  environments 
are  cluttered  and,  therefore,  heterogeneous. 

The  main  focus  of  this  project  was  the  so-called  convexification  method.  This  is  the  globally  convergent  algorithm  of  the  first  generation. 
This  method  was  hilly  investigated.  The  first  breakthrough  result  on  the  convexification  was  reported  in  the  Annual  report  of  2006  and  was 
published  in  2007.  This  publication  got  quite  a  warm  reception  of  the  scientific  community.  Because  the  convexification  is  a  new  method,  it 
is  natural  that  a  number  of  its  different  aspects  was  studied,  which  was  done  in  this  project.  Three  versions  of  the  numerical  realization  of 
the  convexification  were  implemented  and  tested.  Applications  to  imaging  of  both  antipersonnel  land  mines  and  targets  on  battlefields 
covered  by  smogs  and  flames  were  addressed.  In  late  2007  the  second  breakthrough  result  was  obtained.  This  is  a  globally  convergent 
numerical  method  of  the  second  generation.  This  technique  deserves  to  he  investigated  further,  because  it  is  very  promising.  This  method  is 
radically  different  from  the  convexification. 

The  PI  believes  that  globally  convergent  algorithms  for  CIPs,  which  are  developed  in  this  project,  have  a  serious  potential  to  radically 
improve}  the  perfonnance  of  many  imaging  modalities  of  the  interest  to  the  Army.  Along  with  numerical  results,  a  number  of  analytical 
results  were  also  obtained  in  this  project.  Nineteen  (19)  papers  in  refereed  journals  with  the  acknowledgment  of  this  grant  support  were 
published/accepted/ submitted.  Results  of  this  effort  were  presented  at  eighteen  (18)  international  professional  meetings.  One  Ph.D.  thesis 
was  defended  based  on  some  results  of  this  project. 
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ABSTRACT 

Coefficient  Inverse  Problems  (CIPs)  for  Partial  Differential  Equations  (PDEs)  represent 
a  very  important  tool  for  such  needs  of  the  Army  as  imaging  of  unknown  targets  hidden 
in  cluttered  heterogeneous  backgrounds.  The  goal  of  this  project  is  the  development  of 
globally  convergent  numerical  methods  for  a  wide  class  of  CIPs.  These  methods  are  tested 
on  mathematical  models  of  the  interest  to  the  Army  such  as  imaging  of  antipersonnel  land 
mines  and  targets  on  battlefields  covered  by  smogs  and  flames.  In  our  definition  “global 
convergence”  entails:  (1)  a  rigorous  convergence  analysis  that  does  not  depend  on  the  quality 
of  the  initial  guess,  and  (2)  numerical  simulations  that  confirm  the  advertised  convergence 
property. 

A  conventional  way  to  solve  a  CIP  is  via  the  minimization  of  a  least  squares  objective 
functional.  This  functional  characterizes  misfit  between  the  data  and  the  solution  of  that  PDE 
with  a  “guess”  for  the  unknown  coefficient.  However,  it  is  well  known  to  researchers  working 
on  computations  of  inverse  problems  that  the  phenomenon  of  multiple  local  minima  of  these 
functionals  represents  the  major  obstacle  for  the  development  of  reliable  numerical  methods 
for  multidimensional  CIPs.  This  phenomenon  in  turn  is  caused  by  the  above  mentioned  non¬ 
linearity  and  ill-posedness.  Because  of  local  minima,  one  should  somehow  guess  in  advance 
about  a  good  approximation  for  the  solution.  Without  the  availability  of  a  first  good  guess, 
however,  there  is  no  guarantee  that  the  calculated  coefficient  is  indeed  close  to  the  correct 
one.  In  our  terminology  these  are  locally  convergent  numerical  methods.  In  other  words,  their 
convergence  to  the  correct  solution  can  be  guaranteed  only  if  the  starting  point  is  located 
in  a  small  neighborhood  of  this  solution.  Because  of  local  minima,  conventional  numerical 
methods  for  multidimensional  CIPs  are  locally  convergent  ones.  However,  in  many  important 
applications  the  first  good  guess  is  unavailable.  In  particular,  locally  convergent  algorithms 
are  fundamentally  unsatisfactory  for  the  needs  of  the  Army,  because  an  accurate  a  priori 
knowledge  of  the  properties  of  a  medium  is  rarely  available  in  military  applications.  This  is 
because  military  environments  are  cluttered  and,  therefore,  heterogeneous. 

The  main  focus  of  this  project  was  the  so-called  convexification  method.  This  is  the 
globally  convergent  algorithm  of  the  first  generation.  This  method  was  fully  investigated. 
The  first  breakthrough  result  on  the  convexification  was  reported  in  the  Annual  report  of 
2006  and  was  published  in  2007.  This  publication  got  quite  a  warm  reception  of  the  scientific 
community  (section  1).  Because  the  convexification  is  a  new  method,  it  is  natural  that  a 
number  of  its  different  aspects  was  studied,  which  was  done  in  this  project.  Three  versions  of 
the  numerical  realization  of  the  convexification  were  implemented  and  tested.  Applications 
to  imaging  of  both  antipersonnel  land  mines  and  targets  on  battlefields  covered  by  smogs 


1 


and  flames  were  addressed.  In  late  2007  the  second  breakthrough  result  was  obtained. 
This  is  a  globally  convergent  numerical  method  of  the  second  generation.  This  technique 
deserves  to  be  investigated  further,  because  it  is  very  promising.  This  method  is  radically 
different  from  the  convexification. 

The  PI  believes  that  globally  convergent  algorithms  for  CIPs,  which  are  developed  in  this 
project,  have  a  serious  potential  to  radically  improve  the  performance  of  many  imaging 
modalities  of  the  interest  to  the  Army.  Along  with  numerical  results,  a  number  of  analytical 
results  were  also  obtained  in  this  project.  Nineteen  (19)  papers  in  refereed  journals  with  the 
acknowledgment  of  this  grant  support  were  published/accepted/submitted.  Results  of  this 
effort  were  presented  at  eighteen  (18)  international  professional  meetings.  One  Ph.D.  thesis 
was  defended  based  on  some  results  of  this  project. 
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1  Introduction 


In  this  final  report  only  most  interesting  results  of  the  project  are  presented.  Other  results 
can  be  found  in  2006  and  2007  annual  reports  as  well  as  in  publications  of  the  PI,  which  are 
cited  in  section  named  “Publications”. 

1.1  State  of  the  art  of  numerical  methods  for  CIPs  and  crucial 
role  of  globally  convergent  numerical  methods 

CIPs  play  a  rapidly  growing  role  in  military  applications.  Indeed,  one  of  the  needs  of  the 
Army  is  to  detect  and  image  unknown  targets.  Examples  of  those  targets  include  land 
mines,  underground  bunkers,  tanks  on  battlefields  covered  by  smogs  and  flames,  etc..  In  all 
these  scenarios  targets  are  incorporated  in  cluttered,  heterogenous  and,  therefore,  unknown 
backgrounds.  Probing  radiations  are  usually  thought  as  electric  and  acoustic  waves  for 
the  first  two  applications  and  light  originated  by  lasers  in  the  third.  Output  radiations 
are  measured  by  sensors.  Interestingly,  the  diffuse  (because  of  smogs  and  flames)  light 
propagation  in  the  third  application  is  even  helpful.  This  is  because  even  if  the  direct  light 
would  miss  the  target,  diffuse  photons  would  not,  and,  therefore,  detectors  would  still  sense 
the  presence  of  that  target.  Propagations  of  those  signals  are  covered  by  Partial  Differential 
Equations  (PDEs),  which  are  derived  from  the  fundamental  laws  of  physics.  Electric,  acoustic 
or  light  scattering  properties  of  both  unknown  targets  and  the  backgrounds  are  described 
by  coefficients  of  those  PDEs.  Since  such  properties  of  targets  usually  differ  sharply  from 
the  properties  of  the  surroundings,  the  presence  of  targets  of  interest  is  revealed  by  the 
differences  in  the  output  signals  measured  by  sensors. 

Therefore,  the  goal  of  a  CIP  is  to  calculate  a  good  approximation  to  the  unknown  co¬ 
efficient  (s)  of  that  PDE  from  the  measured  data.  However,  the  latter  is  an  enormously 
challenging  mathematical  problem.  Two  main  factors  causing  a  substantial  difficulty  of 
construction  of  stable  globally  convergent  algorithms  for  CIPs  are  their  non-linearity  and 
ill-posedness.  The  nonlinearity  is  because  solutions  of  PDEs  depend  nonlinearly  on  their 
coefficients.  The  ill-posedness  is  a  well  known  feature  of  inverse  problems.  This  means  that 
small  fluctuations  in  the  input  data  can  cause  large  fluctuations  of  solutions.  This  feature 
led  to  the  development  of  regularization  algorithms,  see,  e.g.,  [50]. 

A  conventional  way  to  solve  a  CIP  is  via  the  minimization  of  a  least  squares  objective 
functional.  This  functional  characterizes  misfit  between  the  data  and  the  solution  of  that  PDE 
with  a  “guess”  for  the  unknown  coefficient.  However,  it  is  well  known  to  researchers  working 
on  computations  of  inverse  problems  that  the  phenomenon  of  multiple  local  minima  of  these 
functionals  represents  the  major  obstacle  for  the  development  of  reliable  numerical  methods 
for  multidimensional  CIPs.  This  phenomenon  in  turn  is  caused  by  the  above  mentioned  non¬ 
linearity  and  ill-posedness.  Because  of  local  minima,  one  should  somehow  guess  in  advance 
about  a  good  approximation  for  the  solution.  Without  the  availability  of  a  first  good  guess, 
however,  there  is  no  guarantee  that  the  calculated  coefficient  is  indeed  close  to  the  correct 
one.  In  our  terminology  these  are  locally  convergent  numerical  methods.  In  other  words,  their 


5 


convergence  to  the  correct  solution  can  be  guaranteed  only  if  the  starting  point  is  located 
in  a  small  neighborhood  of  this  solution.  Because  of  local  minima,  conventional  numerical 
methods  for  multidimensional  CIPs  are  locally  convergent  ones.  However,  in  many  important 
applications  the  first  good  guess  is  unavailable.  In  particular,  locally  convergent  algorithms 
are  fundamentally  unsatisfactory  for  the  needs  of  the  Army,  because  an  accurate  a  priori 
knowledge  of  the  properties  of  a  medium  is  rarely  available  in  military  applications.  This  is 
because  military  environments  are  cluttered  and,  therefore,  heterogeneous. 

In  our  definition  “global  convergence”  entails:  (1)  a  rigorous  convergence  analysis  that 
does  not  depend  on  the  quality  of  the  initial  guess,  and  (2)  numerical  simulations  that  confirm 
the  advertised  convergence  property.  Two  new  globally  convergent  algorithms  for  a  broad 
class  of  CIPs  were  developed  in  this  project:  the  convexihcation  [1-7]  and  an  algorithm  which 
is  based  on  a  layer  stripping  procedure  with  respect  to  the  so-called  pseudo  frequency  [12], 
which  is  the  positive  parameter  s  of  the  Laplace  transform  of  either  hyperbolic  or  parabolic 
PDE.  We  call  s  pseudo  frequency.  These  algorithms  represent  respectively  the  first  and  the 
second  generation  of  globally  convergent  numerical  methods.  While  the  convexihcation  was 
fully  investigated  in  this  project,  the  second  one  is  only  in  an  infant  age  and  requires  a 
detailed  further  investigation. 

Since  the  convexihcation  is  a  new  method,  it  is  natural  to  investigate  its  main  features 
from  a  number  of  different  perspectives,  which  was  the  main  goal  of  this  project  and  was 
rehected  in  publications  [1-7].  However,  during  the  work  on  this  project  a  new  main  goal 
occurred  in  2007  [12].  This  one  was  the  development  of  a  new  globally  convergent  algorithm 
for  CIPs,  which  is  based  on  the  layer  stripping  procedure  with  respect  to  the  pseudo  frequency 
rather  than  the  layer  stripping  with  respect  to  a  spatial  variable  (convexihcation).  Since  the 
differential  operator  with  respect  to  the  pseudo  frequency  is  not  a  part  of  a  corresponding 
PDE,  it  is  anticipated  that  this  new  method  will  have  a  good  stability  property.  Thus,  the 
recent  work  [12]  has  started  the  second  generation  of  globally  convergent  numerical  methods 
for  CIPs.  In  addition  to  computational  results,  a  number  of  analytical  results  concerning 
with  the  important  issues  of  uniqueness  and  stability  of  inverse  problems  were  obtained 
during  the  work  on  this  grant  [14-19]. 

1.2  Warm  reception  by  the  scientific  community 

Two  publications  with  the  acknowledgment  of  the  support  of  this  grant  were  publicly  warmly 
accepted  by  the  scientific  community.  First,  this  was  the  result  of  [1],  where  the  first  nu¬ 
merical  result  of  the  convexihcation  in  the  2-d  case  was  published.  As  a  clear  manifestation 
of  a  warm  reception  of  the  publication  [1],  the  PI  was  requested  for  an  interview  by  The 
Institute  of  Physics  Publishing,  www.iop.org,  The  Publisher,  which  publishes  the  journal 
named  “Inverse  Problems”,  which  is  the  most  popular  journal  in  the  held.  Pi’s  interview 
can  be  found  at  http://www.iop.org/ej/authors_edition/.  It  is  stated  at  that  site  that  u60 
Seconds  With  showcases  interviews  with  IOP  authors  who  have  published  papers  that  were 
key  to  the  advancement  of  physics  research  in  their  particular  subject  area .”  In  addition,  a 
theoretical  result  of  the  PI  [18],  which  is  reported  in  subsection  11.2,  was  highlighted  by  the 
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Editorial  Board  of  Inverse  Problems  among  best  publications  of  2006. 


2  An  Outline  of  the  Convexification 

The  convexification  works  with  a  single  location  of  the  source  (equivalently,  a  single  direction 
of  the  incident  plane  wave)  and  with  the  data  collected  at  a  piece  of  the  boundary  rather 
than  on  the  whole  boundary.  This  means  the  minimal  requirement  on  the  data  collection 
scheme.  The  so-called  “transformation  procedure”  of  the  convexification  has  deep  roots  in 
the  method  of  Carleman  estimates  for  CIPs,  which  was  originally  introduced  in  1981  in 
the  work  of  Bukhgeim  and  Klibanov  [29]  (also,  see,  e.g.,  [38],  [39],  [41]  for  some  follow  up 
works)  and  became  since  then  one  of  very  few  classic  tools  in  the  field  of  Inverse  Problems. 
Initially  Carleman  estimates  were  introduced  in  the  field  of  Inverse  Problems  only  for  the 
proofs  of  uniqueness  and  stability  results.  Hence,  prior  to  the  convexification,  Carleman 
estimates  were  not  applied  to  numerical  developments.  The  sequence  of  Carleman  Weight 
Functions  (CWFs)  exp  [—A  (z  —  Zj_i)] ,  i  —  1,  •  •  •  ,  n  for  the  operator  d2/dz2  is  involved  in 
the  numerical  scheme  of  the  convexification.  It  was  explained  in  [1]  that  this  sequence  of 
weights  ensures  the  strict  convexity  of  the  sequence  of  residual  least  squares  functionals,  as 
well  as  the  stabilization  of  the  resulting  layer  stripping  procedure,  also  see  subsection  2.6. 
Here  n  is  the  number  of  layers  {z  :  z  G  (zj_ i,  z,]  } . 

First,  the  original  CIP  is  approximated  with  the  Cauchy  problem  for  a  coupled  system  of 
ordinary  nonlinear  integral-differential  equations  for  the  vector- valued  function  p(z,s).  An 
important  feature  of  this  system  is  that  the  unknown  coefficient  is  not  present  in  it.  From 
the  computational  standpoint  the  major  complicating  factor  of  this  system  is  the  presence 
of  nonlinearities  with  Volterra-like  integrals  with  respect  to  the  pseudo  frequency  s,  whereas 
derivatives  are  taken  with  respect  to  the  spatial  variable  z.  To  solve  this  Cauchy  problem, 
a  stable  layer  stripping  procedure  is  applied.  On  each  thin  z-layer  one  approximates  the 
vector-valued  function  p(z,  s)  as  a  quadratic  polynomial  with  the  unknown  quadratic  term 
and  two  other  terms  being  known  from  the  data  at  the  interface.  Next,  one  constructs 
a  residual  least  squares  functional  J\,i  with  the  CWF  in  it.  The  CWF  for  any  differential 
operator  ensures  that  in  the  weighted  L2  norm  of  this  operator  the  principal  part  dominates 
the  rest.  Hence,  due  to  the  presence  of  the  CWF  in  the  functional  J\.i,  the  principal  linear 
part  p"  of  the  differential  operator  for  the  vector-valued  function  p  in  Ja,;  dominates  all  other 
terms,  including  nonlinear  ones.  Because  of  this  domination,  the  functional  J\tl  is  strictly 
convex  on  a  certain  a  priori  given  bounded  set  (this  is  not  a  small  set!).  Furthermore,  the 
unique  minimizer  of  J\^  belongs  to  the  interior  of  this  set.  In  addition,  this  minimize!'  is 
close  to  the  one  which  corresponds  to  the  exact  solution  of  the  CIP  (Theorem  2.1). 

2.1  Forward  problems 

Below  in  section  2  x  =  (xi,x2,z).  First,  consider  the  Cauchy  problem  for  a  hyperbolic 
equation 

c  (. x )  utt  =  A u  —  a{x)u  in  M3  x  (0,  oo) ,  (2.1) 
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(2.2) 


u  (x,  0)  =  0,  ut  (x,  0 )  —  S  (x  —  xo)  • 

In  addition  we  consider  the  Cauchy  problem  for  a  parabolic  equation 

c  (x)  ut  =  A u  —  a(x)u  in  M3  x  (0,  oo) , 
u  (x,  0)  =  5  (x  —  x0 ) . 

Consider  the  Laplace  transform  of  both  functions  u  and  u, 


w(x ,  s) 


OO 

J  u(x.t)e~stdt 
0 


oo 

J  u(x,t)e~*2tdt 
0 


(2.3) 

(2.4) 


(2.5) 


for  positive  s  >  so  >  0.  Since  both  integrals  of  (2.5)  lead  to  the  same  elliptic  equation  (2.9a) 
for  the  function  w,  it  is  more  convenient  to  us  to  consider  the  theory  for  the  function  w. 

We  use  the  Laplace  transform  with  the  positive  parameter  s  only  because  we  need  to 
make  sure  that  the  function  w(x,s )  >  0  by  the  maximum  principle,  see  (2.10).  As  to  the 
coefficients  of  equations  (2.1)  and  (2.3),  we  assume  that 


c(x)  G  [di,2d2] ,  where  di,d,2  =  const.  >  0.  (2.6) 

c  ( x )  G  C2  (M3)  ,c(x)  =  2di  for  x  G  M3\f2,  (2.7) 

a(x)  G  C 2  (M3)  ,  a(x)  >  0  and  a(x)  =  0  for  x  G  M3\fl,  (2.8) 

where  C  M3  is  a  convex  bounded  domain.  Note  that  in  the  held  of  CIPs  a  certain  over¬ 
smoothness  of  coefficients  is  usually  assumed. 

The  equation  for  the  function  w  is 


Aw  —  [sV2  (x)  +  a(x)]  w  =  —5  (x  —  xo) ,  Vs  >  so  =  const.  >  0.  (2.9a) 

with  the  condition  at  the  infinity 

lim  w(x,  s )  =  0,  Vs  >  s0  =  const.  >  0.  (2.9b) 

|  x\  — >oo 

The  maximum  principle  and  conditions  (2.6)-(2.8)  imply  that  for  all  s  >  So  there  exists 
unique  solution  w(,  s)  G  C3  (M3\  {|a;  —  x0|  <  e}) ,  We  >  0  of  the  problem  (2.9a,b).  Further¬ 
more, 

w(x,s )  >  0,  Vs  >  so-  (2.10) 

Under  certain  conditions  imposed  on  coefficients  c(x)  and  a(x)  the  following  asymptotic 
behavior  holds  [41],  [44] 


1  +  0  (  - 
s 


,  s  — >  oo,  |o;|  <  2,  (3  —  0, 1. 


(2.11) 


This  asymptotic  behavior  was  derived  in  [41],  [44]  on  the  basis  of  Theorem  4.1  in  [47]. 


2.2  The  inverse  problem 

We  formulate  the  inverse  problem  for  the  elliptic  equation  (2.9a)  with  the  condition  (2.9b), 
because  formulations  of  this  problem  for  above  hyperbolic  and  parabolic  equations  are  similar 
due  to  (2.5).  In  the  case  when  the  data  either  for  the  inverse  parabolic  or  for  the  inverse 
hyperbolic  problem  are  given  on  a  large  but  finite  time  interval  t  G  (0,  T)  one  should  assume 
that  in  (2.5) 

OO  T  OO  T 

u(x,t)e~s2tdt,  J  u(x,t)e~stdt  &  j  u(x,t)e~stdt,  s  >  s0  =  const.  >  0. 
o  ooo 

Inverse  Problem.  Let  fl  be  a  rectangular  prism, 

=  {—A  <  xi,  X2  <  A,  z  E  (0,  L )}  , 

where  A  and  L  are  positive  numbers.  Suppose  that  one  of  coefficients  of  the  equation  (2.9a) 
is  unknown  in  and  known  in  M3\f2  and  all  other  coefficients  are  known  everywhere. 
Determine  that  unknown  coefficient  for  G  D,  assuming  that  the  following  two  functions 
if  (aq,  X2,  s)  and  ip  (aq,  X2,  s)  are  known  for  a  single  source  position  x0  ^ 

w  (aq,  x2 ,  0,  s)  =  <f  (aq,  x2,  s ) ,  wz  (x1,  x2,  0,  s)  =  0  (xi,  x2,  s) ,  (2.12) 

for  (xux2 ,s)  G  {-A,  A)2  x  [s0,s], 

where  So  and  s  are  certain  positive  numbers. 

In  the  case  when  only  one  of  functions  <p  or  0  is  given  for  all  (aq ,x2,s)  G  M2  x  (so,  oo), 
another  one  can  be  determined  uniquely  via  solution  of  the  corresponding  boundary  value 
problem  for  the  equation  (2.9a)  in  the  lower  half  space  {z  <  0}  . 

2.3  Transformation 

Since  we  present  numerical  examples  for  the  case  a(x)  =  0,  we  consider  this  case  only  below 
for  brevity.  The  numerical  scheme  for  the  case  of  the  unknown  function  a  (x)  was  considered 
in  [41]  and  [44],  We  first  transform  our  problem  to  the  Cauchy  problem  for  a  nonlinear  elliptic 
integral-differential  equation,  in  which  the  unknown  coefficient  is  not  present.  Because  of 
(2.10)  we  can  consider  the  function  v  =  In  w.  Then  (2.9a)  and  (2.12)  lead  to 

Av  +  |Vu|2  =  s2c(x )  in  D,  (2-13) 

v  (aq,  x2,  0,  s)  =  </q  (aq,  aq,  s) ,  vz  (x1,x2,  0,  s)  =  i/q  (aq,  x2,  s) ,  (2.14) 

for  (aq,  aq,  s)  G  (-A,  A)2  x  (s0,  oo) , 

ip 

<Pi  =  =  -• 

<P 


where 
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Divide  both  sides  of  (2.13)  by  s2  and  denote 


Then 


v(x,  s) 


It  follows  from  (2.11)  that 


Av  +  s 2  |Vx|2  =  c(x) . 


D%v{x,  s)  —  O 


,  D“dsv{x,  s ) 


for 


(X). 


(2.15) 


(2.16) 


We  now  want  to  eliminate  the  unknown  coefficient  c(x)  from  equation  (2.15).  To  do  so, 
differentiate  its  both  sides  with  respect  to  the  parameter  s  and  observe  that  dsc(x)  =  0. 
Denote 


q  (x,  s )  =  dsv  (x,  s ) ,  (p2  (x,  s )  =  ds  [s  2ipl  (x,  s)]  ,  ^2(x^  s )  =  [s  2i>i  (x,  s)]  .  (2.17) 


Then  by  (2.11)  and  (2.16) 

q  (x,  r)  dr. 

We  truncate  this  integral  as 


v  {x.  s)  a;  - — /  q  (x,  t)  dr  +  V(x),  (2-18) 

S 

where  s  >  So  is  a  large  number  which  should  be  chosen  in  numerical  experiments.  This 
truncation  is  similar  to  the  truncation  of  high  frequencies,  which  is  routinely  done  in  science 
and  engineering.  Here  V (x)  is  the  so-called  “tail  function”  which  complements  the  integral. 
This  function  is  unknown  and  V (x)  ~  v  (x,  s) .  Thus,  we  obtain  the  following  (approximate) 
nonlinear  integral  differential  equation 


A  q  —  2  s2Vg 


s 

/» 

S 

/* 

/  \7q  (x,  r)  dr  +  2s 

/  Wq  (x,  r)  dr 

J 

i/ 

S 

_S 

S 


S 


In  addition,  the  following  Cauchy  data  are  given 

q  (xi,  x2,  0,  s)  =  tp2  (xi,  x2,  s ) ,  qz  (xi,  x2,  0,  s)  =  (xi,  x2,  s) , 


(2.20a) 


(2.20b) 
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for  (xi,x2,s)  G  ( -4,4 )2  x  (s0,s). 

Equation  (2.20a)  has  two  unknown  functions  in  it  q  and  V .  A  good  point  of  this  equation, 
however  is  that  it  does  not  contain  the  unknown  coefficient.  Thus,  these  two  unknown 
functions  are  approximated  differently  in  our  algorithms.  One  of  the  key  points  which  helps 
us  to  approximate  the  tail  function  is  that  it  is  small  for  large  values  of  s,  as  it  follows 
from  (2.16).  If  the  tail  function  is  given,  then  the  problem  (2.20a,b)  is  the  Cauchy  problem 
for  a  nonlinear  integral  differential  equation  with  Volterra-like  integrals  depending  on  the 
parameter  s,  which  is  not  involved  in  the  differential  operator.  If  the  second  and  third  terms 
in  (2.20a)  would  be  absent,  then  this  would  be  a  well-known  Cauchy  problem  for  the  Laplace 
equation,  which  is  discussed  in  many  publications  (we  are  not  in  a  position  to  list  those). 
Compared  with  the  latter,  two  major  difficulties  of  the  problem  (2.20a,b)  are  the  nonlinearity 
and  the  presence  of  integrals. 

2.4  Approximation 

Now  the  main  question  is:  How  to  solve  numerically  the  problem  ( 2.20a,b )?  This  question 
is  addressed  in  follow  up  subsections  of  this  section.  Let  {(f) j  (xi,  C  C 2  ([—4,  A]2)  be 

a  set  of  linearly  independent  functions  approximating  the  function  q  (x,  s )  up  to  its  second 
derivatives 

N 

D*q(x1,x2,z,s)  «  D^^Pj  (z,s)(f)j  (x1,x2) ,  (x,  s)  G  Q  x  [s0,  s] ,  \a\  <  2,  (2.21) 

3= i 

We  assume  that  functions  0  -  (xi,x2)  are  given  via  analytical  formulas,  so  that  their  deriva¬ 
tives  can  be  calculated  analytically.  For  example,  one  can  take  functions  (f>]  in  the  form 
4>j  (x\,x2)  =  Qj1  (xi)  Qj2  (. x2 ) ,  where  QJt  and  Qn  are  either  parts  of  an  orthogonal  basis  in 
L2  or  cubic  -B-splines.  Denote  for  brevity  y  :=  (xi,  x2) ,  0  =  (—4,  4)  x  (—4, 4) .  Substitute 
(2.21)  in  (2.20a).  Denote  p'j  ( z,s )  :=  dzpj  ( z,s ) .  Next,  multiply  both  sides  of  the  resulting 
equation  by  the  function  (f>r  (y) ,  r  =  l,---  ,N,  and  integrate  over  the  rectangle  0.  We 
obtain  the  following  coupled  system  of  nonlinear  ordinary  integral-differential  equations 

p  (z,  t)  dr,  =0,  (2.22) 

(2,s)  e  [0,  L]  x  [s0,s] , 

where  B  is  invertible  matrix  and  the  W-dimensional  vector-valued  function  Fed2  (M5iV+1) . 
In  addition  (2.20b)  and  (2.21)  imply  that  the  following  Cauchy  data  p°(s)  and  p1(s) 

P( 0,  s)  =  p°(s),  p'( 0,  s)  =  p\s ),  s  G  [s0,  s] .  (2.23) 

The  main  question  now  is:  How  to  solve  the  problem  (2.22),  (2.23)  in  a  stable  way?  A 
stable  layer  stripping  procedure  is  described  in  the  next  two  subsections. 


II  (p,  z,  s)  Bp"  —  F  p.p,l  p  (z.  t)  dr , 
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2.5  Numerical  Method  for  the  Problem  (2.22),  (2.23) 

Consider  a  uniform  partition  of  the  interval  [0,  L\  with  the  grid  step  size  h, 

0  =  z0  <  Z!  <  ■  ■  ■  <  zn  =  L,  h  =  Zi-  Zi _i. 

In  each  layer  (zi-i,Zi]  we  approximate  the  vector-valued  function  p(z,  s)  with  quadratic 
polynomials 

p(z,  s)  pi(z,  s)  =  — — cii(s )  +  (z-  z^ i)  bi(s)  +  q(s),  (2.24) 

Z  e  ,  s  G  [s0,s] . 

Hence,  q(s)  =  pj(zj_i,  s)  and  6j(s)  =  s).  By  (2.23)  functions  &o(s)  =  Pi(^ 'o,  s)  =  p1  (s) 

and  Co(s)  =  Pi(^0)s)  =  P°  (s)  are  known  from  the  available  data  for  the  inverse  problem. 
Functions  6j(s)  and  q(s)  are  assumed  to  be  known  from  the  previous  step  of  the  layer 
stripping  procedure.  Hence,  the  only  unknown  function  in  (2.24)  is  the  quadratic  term  aj(s). 
As  soon  as  this  term  is  approximately  found,  one  sets 

h2 

Pi{zi,s)  :=  Cj+i(s)  =  yat:(s)  +  hbi(s)  +  q(s), 

Pi(zi,s )  :=  bi+1(s)  =  ha.i(s)  +  bi(s). 

We  now  focus  on  the  procedure  of  hireling  an  approximation  for  the  quadratic  term  oq(s). 
Consider  the  weight  function  C\:i  (z) , 


Cx,i  (z)  =  exp  [—A  (z  -  Zi-i)\ , 


where  A  >  1  is  a  parameter  which  will  be  chosen  later.  It  was  proven  in  [41]  that  C\^  (z)  is  the 
weight  function  involved  in  the  Carleman  estimate  for  the  operator  d2/dz2  on  the  interval 
(zi-i,Zi).  In  other  words,  this  is  a  CWF.  Construct  the  weighted  least  squares  objective 
function  J\ti, 


J\,i  (ai(s))  —  ds  n2  \pi(z,s),s]Cx,i(z)  dz. 


(2.25) 


so 


Zi- 1 


Let  M  be  a  positive  number.  Consider  the  bounded  set  G  (. M )  C  C  [s0,  s] , 

G  (M)  =  |a(s)  G  C[s0,  s]  :  lla(s)llc[s0,s]  — 


(2.26) 


We  search  for  a  minimizer  di(s)  G  G  (M)  of  the  functional  Hence,  (2.24)  and  (2.26) 
imply  that  j^k\z,  s )  <  M',  k  —  0,1,  2,  where  M'  =  M'  (M)  is  a  positive  constant  depending 

on  M.  Hence,  by  the  Archcla  theorem  the  condition  (2.26)  implies  that  the  function  pt(z,  s) 
belongs  to  a  subset  of  a  priori  chosen  compact  set  (depending  on  M)  in  C1  [^-i,  zi\  x 
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C  [so,  s]  ,  which  correspond  well  with  the  above  mentioned  Tikhonov  theorem.  Note  that  our 
theory  does  not  impose  a  “smallness”  condition  on  the  constant  M.  In  accordance  with  the 
Tikhonov  concept  the  constant  M  should  be  chosen  for  a  specific  set  of  applied  problems 
under  consideration.  It  was  proven  in  [41]- [44]  that  the  functional  J\: j  is  strictly  convex 
on  the  set  G  ( M )  for  an  appropriate  choice  of  a  sufficiently  large  parameter  A  >  Ao  (M) . 
Furthermore,  its  unique  minimizer  can  be  found  as  the  unique  solution  of  an  equation  with 
a  contractual  mapping  operator.  In  addition,  the  convergence  theorem  of  [41]- [44]  implies 
that  this  minimizer  is  close  to  the  value  a*  (s)  =  (p*)  ( Zi,s ) ,  where  the  function  p*  (z,s) 
corresponds  to  the  exact  solution.  We  combine  these  statements  in  Theorem  2.1  below. 

We  now  formulate  convergence  theorem  for  our  method  assuming  that  tail  functions 
are  small.  Suppose  that  there  exists  the  exact  solution  q*  (x,  s )  G  C3  (fi  x  [s0,  s])  of  the 
equation  (2.20a)  with  the  Cauchy  data  (2.20b).  We  assume  first  that  these  Cauchy  data 
are  “ideal”  ones,  i.e.,  they  are  given  without  an  error.  If  such  a  solution  exists,  then  it  is 
unique,  see  Theorem  6.5.1  in  [41].  Let  p*  (z,  s )  be  the  vector- valued  function  p  (z,  s )  which  is 
obtained  from  q*  (x,s)  via  the  approximation  (2.22).  Since  (2.22)  is  only  an  approximation, 
we  assume  that  the  function  p*  ( z,s )  G  C3([0,L]  x  [so,s])  satisfies  the  following  analog  of 
equation  (2.24)  (recall  that  we  set  tails  to  zero  in  our  convergence  theorem) 


B  {p*)"  (z,  s)  —  F  (p*)'  (z,s)  ,p*(z,s),  /  (p*)'  (z,r)  dr,  /  p*  (z,  r)  dr,  0,  0,  s  (2.27) 


=  e(z,s),  z  G  (0,  L) ,  s  G  [s0,  s] . 
And  by  (2.23)  Cauchy  data  for  p*  (z,  s )  are 

p*  (0,  s)  =  p°*  (s) ,  (p*)'  (0,  s)  =  p1*  (s) . 


(2.28) 


We  assume  that  the  function  £  (z,  s )  is  sufficiently  small  (actually,  this  function  is  unknown), 


lk(^»s)llc([o,L]x[ao,5])  -£’  (2-29) 

where  £  is  a  small  positive  number.  Because  in  the  reality  the  Cauchy  data  (2.23)  are  always 
given  with  an  error,  we  assume  that 

V  (S)  -  P° *  (S)IUo,s]  +  Wpl  ^-P1*  (S)|lcp0.s1  ^  (2-3°) 

Denote 

h 

1  _  p-A/i  r 

Iq  (A,  h)  =  - - - =  /  e~Xzdz.  (2.31) 

o 

Let  G  (. M )  be  the  set  of  functions  a  (s)  G  G  (. M )  with  the  L2  (so,  s )  —norm  in  it.  Hence,  by 
(2.26) 

aeC  [s0,  s] ,  ||  o  ||  (7(50,8]  <  ll«IL2(S0,s)  ^  MVs  -  s0,  Va  G  G  (M) . 
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We  formulate  Theorem  2.1  for  the  case  of  a  general  coupled  system  of  ordinary  integral 
differential  equations. 

Theorem  2.1  [41]- [44].  Let  p(z,s )  G  C 2  [0,  L\  x  C  [.s-0,  s]  be  an  N -dimensional  vector- 
valued,  function  and  F  G  C2  (M5JV+1)  also  be  an  N -dimensional  vector-valued  function.  Let 
the  vector-valued  function  p*  ( z ,  s )  G  C 3  ([0,  L\  x  C  [s0,  s])  satisfies  conditions  (2.27)-  (2.30) 
and  ||p*  (z,  s)||C3j0  LixCrSo  ^  <  M /2.  Then  there  exists  a  positive  constant  C  =  C  ( M ,  F,  s0,  s,  L ,  N ) 
and  small  positive  numbers  £o  =  £o  (M,  F,  s0,  s,  L,  N) ,  ho  —  ho  ( M ,  F,  s0,  s,  L,  N )  depending 
only  on  numbers  sq,s,L,M  and  the  vector-valued  function  F  such  that  if  £  G  (0,  £o) ,  h  G 
(0,  ho), 

C 

A  >  A0  := 

£ 

and  that  the  tail  function  is  such  that  ||VV||c/^\  <  £,  then  functions 

Pi  (Zi- 1,  s) ,  %  (Zi- 1,  s)gG  (M) ,  i  =  1,  •  •  •  ,  n, 


and  all  functionals  J\^  are  strictly  convex  on  the  set  G  (. M )  with  the  property 


Iq  (A,  h) 


[ J\,i  (a  +  b)  -  Jx,i  (a)  -  J'X  l  (a)  ( b )]  >  p 


1 2 

I  Z/2(«0,s) 


Vo,  a  +  b  e  G  (M) ,  (2.32) 


where  J'x  t  (a)  is  the  Frechet  derivative  of  the  functional  JXti  at  the  point  a,  and  the  convexity 
constant  p  G  (1/2,1).  The  unique  minimizer  a x^  G  G  (M)  of  the  functional  JXti  is  an 
interior  point  of  the  set  G  (. M )  and  can  be  found  as  the  solution  of  the  equation 


a  =  <hA ,i  (a) ,  i  =  1,  •  •  •  ,n 


with  the  contraction  mapping  operator  $^,1  :  G  (M)  — >  G  (. M )  and 

C 

||aAl,i  —  0'A2,i||  <  VAi,  A2  >  Ao,  i  —  1,  •  •  •  ,  n. 
Furthermore,  the  following  stability  and  convergence  estimate  holds 


(2.33) 


(2.34) 


max  sup  \\Drz  [pt  (z,  s )  -  p*  (z,  s)]  ||c([  5])  <  C  (e  +  h) ,  r  =  0, 1,  2.  (2.35) 

Remark.  We  assume  that  the  tail  function  is  sufficiently  small  because  it  is  indeed 
small  for  sufficiently  large  truncation  pseudo  frequency  s,  see  (2.16).  The  function  J0  (A ,h) 
(see  (2.31))  plays  the  role  of  a  calibration  factor  here.  The  estimate  (2.34)  means  that 
the  minimizers  are  “almost”  independent  on  the  parameter  A,  as  long  as  this  parameter  is 
sufficiently  large.  The  estimate  (2.35)  implies  that  the  difference  between  the  exact  solution 
and  the  one  obtained  by  the  above  method  is  small  as  soon  as  parameters  £  and  h  are 
sufficiently  small.  It  is  important  that  estimates  (2.34),  (2.35)  are  actually  confirmed  in  our 
numerical  experiments,  see  section  5. 
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A  peculiar  question  of  Theorem  2.1  is  that  the  strict  convexity  of  functionals  J\y *  is 
guaranteed  only  on  the  bounded  set  G  ( M )  rather  than  on  the  whole  space  L2  (s0,  s) .  Hence, 
it  seems  to  be,  at  least  at  the  first  glance,  that  in  a  practical  use  of  a  gradient-like  method 
one  should  make  sure  that  points  resulting  from  iterations  of  such  a  method  belong  to  the 
interior  of  G  (M) ,  which  would  complicate  things.  However,  Theorems  1.2  and  1.3  of  [40] 
ensure  that  one  should  not  be  concerned  with  the  latter.  These  theorems  basically  state  that 
if  the  exact  solution  belongs  to  the  set  G  (M/2)  and  the  starting  point  is  also  chosen  on  this 
set,  then  the  distance  between  the  sequential  points  obtained  by  the  gradient  method  and 
the  true  minimizer  of  a  strictly  convex  functional  is  decreasing,  np  to  a  certain  level,  as  the 
number  of  iterations  of  the  gradient  method  is  increasing.  This  level,  in  turn  is  determined 
by  the  level  £  of  the  error  in  the  data.  Thus,  assuming  that  our  exact  solution  p*  is  such 
that  (Theorem  2.1)  ||p*  (z,  -s)||c.3q0  nxrso  5j)  <  M/2  and  starting  iterations  of  a  gradient-like 

method  from  a  point  at(s)  G  G  (M/2)  ,  one  is  guaranteed  that  consecutive  iterations  will  not 
lead  outside  of  G  (M). 


2.6  A  crucial  role  of  Carleman  Weight  Functions 

It  is  important  for  the  understanding  of  the  convexification  to  explain  the  role  of  the  Carle- 
man  Weight  Functions.  Let  the  vector  function  p  (z,  s )  G  C3  [0,  L\  x  C  [so,  s\  be  the  solution 
of  the  Cauchy  problem  (2.22),  (2.23).  By  the  Taylor  formula 


p(z,s) 


p"(zi- i,s) 


+  p'  (A-i,s)  (z  -  Zi_ i) 


+  p(zi_1,s) 


+  o  /  p"  (€>  s )  (z  -  O2  dZ,  (a  s)  G  (Zi- 1,  Zi]  x  [s0,  s  . 


Zi- 1 


Hence,  the  quadratic  approximation  (2.24)  provides  O  (z  —  Zi- 1)  error  for  the  equation  (2.22) 
for  (z,s)  G  (zj_i,Zj]  x  [s0,s]  as  (z  —  Zi- 1)  — >  0+.  Recall  that  Zi  —  Zi- \  =  h.  Hence,  if  the 
CWFs  would  not  be  present,  then  after  considering  n  =  L/h  layers,  one  would  have  the  total 
error  of  at  least  0(nh )  =  O(L),  which  is  large.  Denote  E^z,  s )  the  error  function  due  to  the 
approximation  (2.24)  on  the  layer  (zj_i,  Zj],  Hence, 


Ei(z,s)  =  (z  —  Zi-i)  Ei(z,s),  where  Ei(z,s)  <  M  for  (z,  s)  G  [zi-i,Zi]  x  [s0,  s]  ■  (2.36) 


Denote 


Fi(s)  =  Q  bi(s),Ci(s),  bi(r)  dr,  Ci(r)  dr,  Zi-i,s  := 


Q  p  (zi_i,s)  ,p(zi-us),  /  p  (zi-i,r)dT,  /  p(zi-1,r)dr,zi_1,s  . 
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Substituting  (2.24)  in  (2.22)  and  using  Taylor  formula,  one  can  represent  the  equation  (2.22) 
in  the  form 


ai(s)  -  Fi(s)  -  (z  -  Zi-i)  Hi  z,s,ai(s),  ai(r)  dr  —  (z  —  Zj_i)  Et(z,  s)  =  0,  (2.37) 


where  H \  is  a  C1— function  of  its  variables.  One  can  easily  prove  that 

h 


z  \f(z)\  e~Xzdz  <  -  [he~Xh  +  /„( A,  h)]  \\f\\c[m ,  V/  G  C  [0,  h] ,  (2.38) 


where  J0(A, h )  was  introduced  in  (2.31). 

Suppose  that  we  choose  such  values  of  the  parameter  A  that  A h  >  1.  Since  h/X  =  (A h)  /A2, 
then 


(Xhj  J  z\f(z)\e  Xzdz<  ^||/||C[(ylpV/eC'[0>/i]>  (2-39) 


Hence,  by  (2.36) 


Also, 


Zi 

1  [  ( 

h(\h)  ]  Z'~l] 

Ei(z,s ) 

2  M 


Zi-1 


Zi 


Io(X,h) 


(z  -  Zi-])2  Ei(z,s)  Cx,i(z)dz< 


2  M2h 
A  ‘ 


(2.40) 


(2.41) 


Zi-1 


Thus,  it  follows  from  (2.39)-(2.41)  that  after  n  =  L/h  layers  the  error  due  to  the  quadratic 
approximation  will  be  of  the  order  of 


Choosing  A  =  A  (h,  e)  such  that 


0IT- 


A  >  h2  ’ 


we  obtain  that,  regardless  on  the  approximation  error  Ei(z,s )  in  (2.37),  the  estimate  (2.35) 
of  the  convergence  rate  still  holds. 


Summary 


It  is  shown  both  in  this  subsection  and  in  Theorem  2.1  that  the  role  of  Carleman  Weight 
Functions  is  threefold.  They  ensure:  (1)  the  stability  of  the  layer  stripping  procedure,  (2) 
the  convergence  of  the  convexihcation  method  and  (3)  the  strict  convexity  of  the  residual 
functionals  J\.y 
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2.7  Inversion 

Suppose  that  functions 

drzp(z,s)  :=  drzpi(z,s),  ( z,s )  £  (2i_i,Zi]  x  (s0,  s) ,  i  =  1,  •  ■  ■  ,n,  r  =  0, 1,2  (2.36) 

are  found  by  the  above  layer  stripping  procedure  and  tails  are  also  approximated.  To  ap¬ 
proximate  the  target  coefficient  c(x),  we  use  the  backwards  calculation.  First,  we  calculate 
the  function  v(xi,X2,z,s)  and  its  corresponding  derivatives  by  (2.21)  and  (2.18).  Finally, 
we  approximate  the  coefficient  c(x)  via  calculating  the  left  hand  side  of  (2.36)  at  s  :=  s0- 
In  principle,  one  can  use  any  value  of  s  in  (2.36).  However,  our  computational  experience 
shows  that  the  best  value  is  the  lowest  one  s  :=  Sq-  An  explanation  of  this  is  that  at  s  So 
the  truncation  at  s  =  s  of  the  integral  (2.18)  causes  the  least  error. 

2.8  The  first  numerical  result  in  2-D 

We  now  present  the  first  2-D  numerical  result,  which  was  obtained  by  the  convexification, 
and  which  has  impressed  the  scientific  community  (subsection  1.2),  when  it  was  published 
in  [1]  (also,  see  [2]).  This  result  was  obtained  by  Dr.  A.  Timonov.  Although  it  seems  at 
the  first  glance  that  this  result  is  restricted  only  to  a  medical  application,  but  actually  since 
we  work  with  the  back-reflected  data  here,  Army  applications  are  quite  feasible.  The  2-D 
inverse  problem  of  the  determination  of  the  unknown  electric  conductivity  coefficient  a  (x,  y ) 
in  the  equation 

uXx  +  Uzz  —  iupcr  (x,  z)u  —  0, 

was  solved.  Here  u  is  the  modulated  frequency  and  p  =  An  ■  10  ‘  H/m  is  the  magnetic 
permeability  of  the  vacuum.  The  problem  of  microwave  imaging  of  the  human  abdomen 
was  considered.  The  realistic  case  of  the  frequency  sweep  uj  £  2n  ■  [0.5, 10]  GHz  was  used. 
Since  the  maximum  sensitivity  depth  of  the  back-reflected  microwave  signal  is  6  centimeters 
(see  reference  [16]  in  the  annual  report  [21]),  the  X-ray  Computed  Tomography  image  of 
the  abdomen  was  re-scaled  to  6  cm  thickness.  In  data  simulations  for  the  forward  problem 
a  realistic  range  of  the  parameter  a(x,z)  £  [0.4, 4.8]  Siemens/meter  was  used.  Namely, 
different  values  of  this  parameter  were  assigned  to  different  parts  of  the  CT  image:  highest 
values  were  assigned  to  the  white  areas  and  lowest  to  the  white  ones.  Next,  those  values 
were  linearly  interpolated  for  the  rest.  Instead  of  the  point  source  the  initializing  plane  wave 
elu}Z  was  used. 

This  plane  wave  has  propagated  downwards  (positive  direction  of  the  z-axis)  and  then 
was  scattered  when  entering  the  heterogeneous  medium.  Measurements  of  the  back-reflected 
wave  were  modeled  on  the  top  side  of  the  image.  Additive  random  noise  of  1%,  5%  and 
10%  was  added  to  the  data.  The  inverse  problem  was  solved  for  15  realizations  of  the  noisy 
data  and  then  averaged  the  resulting  functions  a  (x,  z) .  Trigonometric  series  were  used  in 
(2.21)  with  (f>j  (. x )  =  sin  (. kaj )  and  also  (pJ  (x)  =  cos  ( kaj ) ,  where  the  parameter  a  depends 
on  the  size  of  the  medium  in  the  horizontal  direction  (to  make  sure  that  these  functions 
are  orthogonal  in  L2).  We  have  used  k  £  [0,  N]  with  N  £  [10,20] .  Significant  differences  in 
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Figure  1:  The  original  (top  left)  and  reconstructed  images  of  the  function  a  ( x ,  z)  via  the 
convexification.  The  plane  wave  falls  from  the  above  and  measurements  of  the  backreflected 
signal  were  also  simulated  at  the  top  side  of  the  square.  The  top  right  image  represents  the 
reconstruction  result  at  the  1%  noise  in  the  data,  the  bottom  left  and  bottom  right  images 
correspond  to  5%  and  10%  noise  respectively. 

images  with  the  change  of  N  were  not  observed.  As  to  the  tail  function  V(z,z),  we  have 
taken  the  one,  which  corresponds  to  the  uniform  background,  i.e.,  the  background  outside 
of  the  heterogeneous  medium.  This  background  was  assumed  to  be  known,  according  to  the 
above  statement  of  the  inverse  problem.  Figures  1  and  2  display  the  computational  results. 

3  The  Second  Numerical  Implementation  of  the  Con¬ 
vexification 

Results  of  this  and  next  sections  represent  the  joint  effort  of  the  Postdoctoral  Research 
Associate  Dr.  J.  Xin  and  the  PI  [3-6].  Dr.  Xin  has  been  working  on  a  new  numerical 
implementation  of  the  convexification.  Dr.  Xin  has  been  fully  supported  by  this  grant  from 
8/15/05  to  6/30/07.  His  effort  has  been  focused  on  the  development  of  new  numerical 
ideas  of  the  convexification  and  their  numerical  implementation.  In  these  publications  a 
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Figure  2:  The  cross-section  of  images  of  Figure  1  by  horizontal  line  through  the  middle. 
Solid  line  corresponds  to  the  original  image.  Dotted  and  dashed  lines  correspond  to  images 
at  1%  and  5%  noise  respectively. 


simplified  mathematical  model  of  imaging  of  the  spatially  distributed  dielectric  constant  in 
antipersonnel  land  mines  was  treated  by  the  convexification. 

Compared  with  previous  works  [1,2]  on  the  convexification  algorithm,  four  new  ingredi¬ 
ents  of  the  second  numerical  implementation  are:  (I)  we  minimize  strictly  convex  functionals 
directly  for  each  generic  layer  rather  than  via  the  solution  of  the  equivalent  equation  (2.33) 
with  the  contraction  mapping  operator,  (II)  a  local  basis  consisting  of  cubic  i?-splincs  is 
employed  in  the  spatial  approximation  instead  of  the  global  trigonometric  basis  [1],  [2], 
hence,  enabling  sharper  resolution  of  the  reconstructed  material  property  at  the  interface 
between  the  inclusion  and  the  background,  (III)  tails  in  truncated  integrals  (2.18)  are  fitted 
in  to  compensate  the  missing  information,  and  (IV)  we  approximate  the  functions  which 
depend  on  the  pseudo-frequency  s  by  Legendre  polynomials,  thus  calculating  the  integrals 
involving  psendo-freqnency  explicitly  rather  than  numerically.  An  advantage  of  the  direct 
minimization  of  the  convex  functionals  over  solving  the  operator  equation  (2.33)  is  that 
the  time-consuming  pre-programming  and  pre-computational  effort  to  derive  the  operator 
equation  is  avoided. 

In  addition  to  the  new  ingredients  of  the  numerical  implementation,  a  systematic,  com¬ 
parative  analysis  of  various  aspects  on  the  convexification  algorithm  was  performed  [4].  No 
such  study  has  been  performed  before,  and  the  parameters  in  the  algorithm  were  taken  in  an 
ad  hoc  manner.  Certainly  the  comparative  study  is  useful  for  the  further  development  of  the 
convexification  algorithm.  Throughout  this  section  we  use  notation  Sz  for  h,  the  thickness 
of  each  layer  in  the  layer  stripping  procedure  of  the  convexification. 
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PERMITTIVITY  AMD  C0MDUCTIV1TY 
(DASHED}  VS.  FREQUENCY  FOR 
SANDY  LOAM  FOR  VARIOUS  MOISTURES 
(BY  PERCENT  DRY  WEIGHT) 
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Figure  3:  Dielectric  permittivity  and  electric  conductivity  constants  in  soil. 

3.1  A  simplified  mathematical  model  of  imaging  of  antipersonnel 
land  mines 

We  present  a  mathematical  model  for  our  Coefficient  Inverse  Problem  -  identification  of 
antipersonnel  land  mines.  We  work  with  a  simplified  2D  model  and  consider  realistic  ranges 
of  parameters.  We  neglect  the  irregularities  of  ground  surface.  In  order  to  avoid  difficulties 
of  solutions  of  forward  and  inverse  problems  at  the  air-ground  interface,  we  assume  that  the 
dielectric  permittivity  £  of  the  media  is  continuous  everywhere.  Further,  we  take  no  account 
of  the  electric  conductivity  of  the  media.  This  can  be  justified  in  the  case  of  dry  soil,  whose 
electric  conductivity  is  small,  see  Figure  3. 

Suppose  a  pulse  generating  a  polarized  electric  field  occurs  at  the  point  Xo(0,  — £)  when 
the  initial  time  is  t  —  0,  where  £  =  const.  >  0.  The  following  hyperbolic  equation  can  be 
derived  from  the  Maxwell  equations 

—fi£(x)utt  +  A u  =  0,  (x,t)  G  M2  x  M+,  (3.1) 
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u  (x,  0)  =  0,  ut  (x,  0 )  —  S  (x  —  xo)  •  (3.2) 

In  (3.1),  (3.2)  the  function  u(x,t)  is  one  component  of  the  electric  field  and  the  parameter 
H  =  4tt  x  10“'  ( He'ilry )  is  the  magnetic  permeability  in  free  space  and 

£  =  £o£r{x)  (3.3) 

is  the  dielectric  permittivity,  where  £o  ~  8.854  x  10“12  (Famd)  is  the  dielectric  permittivity 
of  free  space  and  er(x)  is  the  dimensionless  relative  dielectric  permittivity  of  the  medium. 

In  both  dry  soil  and  trinitrotoluene  (TNT)  we  have  £r  ~  2.9,  see  Figure  3.  We  are 
interested  in  the  identification  of  antipersonnel  plastic  mines,  which  is  more  difficult  in  a 
practical  scenario  since  the  metal  component  in  them  is  very  small.  Thus,  we  need  to  find 
one  parameter  inside  the  mine  which  can  give  us  sufficient  contrast  with  the  surrounding  dry 
soil.  It  is  well  known  that  a  noticeable  part  of  the  volume  of  any  mine  is  filled  with  air  and 
£r  —  1  for  the  air.  Since  the  mine  does  not  wholly  consist  of  air,  it  is  reasonable  to  assume 
that  £r  =  1.5  inside  the  mine,  which  is  about  the  average  value  of  the  coefficient  £r  within 
the  mine.  Therefore,  for  all  our  computation  we  assume 


J  2.9  outside  mines 
1.5  inside  mines 


(3.4) 


The  sizes  of  antipersonnel  mines  vary  between  5  cm  and  10  cm,  and  usually  they  lay  at  a 
small  depth  underneath  the  ground,  not  exceeding  10  cm.  Hence,  we  model  mines  as  disks 
of  radius  5  cm  which  are  located  in  the  range  z  G  [0,  9]  cm. 

If  we  could  identify  the  coefficient  £r(x),  then  points  whose  values  are  close  to  1.5  will 
be  those  inside  or  close  to  the  mine.  Thus,  finding  an  approximation  for  this  coefficient  via 
solution  of  the  inverse  problem  formulated  below  could  provide  one  with  a  useful  information 
about  the  location  and  relative  dielectric  permittivity  of  the  mine.  This  value,  in  turn,  might 
in  the  future  serve  as  one  parameter  in  a  ‘classifier’  procedure,  which  would  distinguish  mines 
from  the  clutter. 

Consider  the  Laplace  transform  of  the  function  w(x,  t) 

OO 

o 

Since  the  parameters  /i  and  £q  are  rather  small,  by  combining  the  parameter  s  we  re-scale 
these  two  parameters  and  introduce  a  new  variable  k  :=  Sy/JH o,  which  we  call  “pseudo¬ 
frequency”.  Equation  (3.1)  with  the  initial  condition  (3.2)  is  transformed  into 

—Aw  +  K2£r(x)w  —  S  (x  —  xo) ,  x  G  M2,  (3.4) 

lim  w(x,  k)  =  0.  (3.5) 

|x|  — »oo 

We  formulate  our  inverse  problem,  termed  IP. 
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Inverse  Problem  (IP).  Consider  a  rectangle  QcK2 

Q  :=  {-A  <  x  <  A,  z  G  (0,  L)} ,  A,L>  0. 

Suppose  the  coefficient  £r(x)  of  the  equation  (3.4)  is  unknown  in  hi  and  is  a  known  positive 
constant  in  M2\f2.  Determine  the  relative  dielectric  permittivity  £r{x)  for  i6fl,  assuming 
the  following  two  functions  ip  (x,  k)  and  ^  (x,  hi)  are  known  for  a  single  source  position  Xq 

w  (. x ,  0 ,k)  —  (p  (. x ,  hi) ,  wz  (x,  0,  k)  —  "0  (x,  hi) ,  V  (x,  hi)  G  {—A,  A)  x  [k0,  R] .  (3.6) 

We  need  to  decide  the  lower  limit  hi0  and  upper  limit  R  of  pseudo-frequency  hi  for  our 
inverse  problem.  To  find  an  appropriate  constant  R,  we  compute  solutions  of  the  forward 
problem  (3.4),  (3.5)  for  different  values  of  the  parameter  n  >  0  and  determine  such  a  value 
k  :=  R,  beyond  which  the  asymptotic  behavior,  i.e. ,  exponential  decay,  of  the  function  w(x,  hi) 
holds,  see  (2.11).  Thus,  we  identify  such  large  values  of  hi  for  which 

lnu;(x,  hi)  =  v(x,  hi)  fa  l1(x)hi  +  l o(x)  (3.7) 

for  many  points  x  G  fl.  The  function  v(x,  hi)  looks  like  a  straight  line  with  respect  to  n  for 
n  >  R.  We  solve  the  forward  problem  on  a  large  domain  S  :=  (—6  <  x,  z  <  6}  using  the 
finite  element  package  -  COMSOL  Multiphysics™  version  3.2.  Zero  Dirichlct  boundary  is 
imposed  on  the  boundary  <9S.  We  use  triangular  elements  with  Lagrange  cubic  basis.  We 
obtain  R  —  10  after  rounding  to  integers.  One  is  free  to  choose  the  lower  bound  Kq .  We  may 
take  a  smaller  value  that  is  close  to  R,  e.g.,  Hq  =  {8,9}.  The  influence  of  the  lower  limit  k0 
on  the  resolution  of  the  reconstructed  unknown  coefficient  £r  was  also  studied  in  [4], 

3.2  Tails 

Recall  that  by  (2.18) 

K, 

v  (x,  n)  =  —  j  q  (x,  r)  dr  +  V  (x)  ,  (3.9) 

where  V  (x)  =  v  (. x ,  R)  is  the  so-called  “tail  function” .  This  function  is  unknown  and  we  now 
describe  the  procedure  of  approximating  it.  While  in  [1,2]  the  tail  function  was  taken  from 
the  uniform  background  (subsection  2.8),  we  use  here  a  different  procedure  to  approximate 
this  function.  First,  we  represent  this  function  similarly  with  (2.21), 

N 

D%V  (xi,  x2:  z)  «  D"  ^  Vj  (z)  <j>j  (xi,x2) ,  iGO,  \a\  <  2.  (3.10) 

j= 1 

Below  in  this  section  we  describe  a  procedure  of  approximating  functions  V3  (z) ,  V-  (z) ,  V"  (z) . 
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3.2.1  Piecewise  constant  approximation 


From  (2.11)  and  (2.15)  the  asymptotic  behavior  of  the  function  v  (x,  hi)  =  In  w  (x,  n)  /  n2  at 
hi  — >  oo  is 

MX)  •  9lW+o( 4),«-oc, 


V  \X.  K)  = 


+ 


hi 


K 


KT 


(3.11) 


where  functions  j\  (x)  and  gi(x)  are  unknown  because  the  coefficient  er  (x)  is  unknown. 
Therefore,  they  should  be  found  approximately.  These  two  functions  are  connected  with 
each  other  by  the  relation  g±(x)  =  —  In  (— 47t/i(x)).  However,  we  do  not  use  this  connection 
explicitly  in  this  study.  It  follows  from  (3.11)  and  (2.21)  that 


—  t  k  — »  oo,  |ok|  <  2,  (3  —  0, 1. 


We  approximate  functions  j\ (x),  g\{x)  sequentially  layer-by-layer  ignoring  the  higher  order 
term  O  (1/k3).  We  approximate  functions  f\  (x),  gi(x)  as  constants  with  respect  to  z  in  each 
thin  layer  z  G  [^_i,2j),  i.e., 

fi(x)  :=  fi(x,y,Zi-i),  gi(x)  :=  gi(x,  y,  ^_i),  z  G  [^_i,^).  (3.12) 

We  also  set 

dzfi(x)  ■=  d%fi(x,y,Zi-i),  d^gi(x)  :=  d%gi(x,y,Zi-i),  z  G  [^_i,^)s  /?  =  1,  2.  (3.13) 

By  (3.9),  (3.10)  and  (2.21) 


N  „  N  N 

V[x,s)  =  -^(pjix^)  /  pj(z,T)dT  +  ^Hj(z)^(x,?/)  =  ^vj(z,s)(f)j(x,y), 
i=i  {  i=i  l=i 


where 

K, 

Vj(z,s)  =  -  J Pj(z,r)dr +  Vj(z).  (3.14) 

Furthermore,  from  (3.11)-(3.13)  we  assume  that  for  a  sufficiently  large  number  k  G  (kq,k) 


Vj  [z,  K)  = 


kj  (Zi—  l)  tj  (-Zj_  l)  r  \  \ 

- 1  ^ — ,  z  G  pi-i,  z?)  ,  kG  (k,  n)  , 


K 


K 


~~  (  \  /C17)  (^-l)  ,  47)  (Z*-l) 

<9>i  (a,  k)  = 


/c 


/c 


,  z  G  [z*-i,  Zj) ,  kG  (k,  k)  ,  7  =  1,2. 


Once  numbers  k ,-7')  (^_i)  are  found,  we  set 


,7  =  0, 1,2. 


(3.15) 

(3.16) 
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Suppose  we  found  all  numbers 

kj  (zi- 1) ,  tj  (zi- 1) ,  kj7)  (zi- 1) ,  tj7)  (zi- 1) . 

Thus,  by  (3.11)-(3.16)  we  set 

N 


D(*,y)v  (*>  K )  :=  X 
3= 1 


-  [ Pj(z,r)dr  +  hhzll  +  ^ 


K 


7c2 


N 


■=  X  vi(z’  (x>  y)  ’  z  G  fc-i»  *0  >  «  e  [«o,  «] ,  M  <  2, 

j=i 


JV 


d%v(x,K)  :=  X 


3=1 


-  /  9]pj(z,r)dr  + 


fcj7)  (fr-i)  +  47)  fe-l) 


AC 


/C2 


AT 


:=  X]  3>j(z,  K)<f>j  (x,  y),  ze  [Zi- 1,  Zi)  ,  AC  G  [/Co,  Ac]  ,  7=1,2. 

J=1 

We  end  up  with  finding  the  approximate  values  in  (3.17). 


(3.17) 


DUM^y)  (3-18) 


(pAx,y)  (3.19) 


3.2.2  Approximation  of  numbers  in  (3.17) 

We  consider  a  step-by-step  procedure  to  find  these  numbers. 

Step  1.  Consider  the  first  layer  z  G  [0,  Zi]  ,  with  z0  —  0,  Z\  —  Sz. 
We  compute  functions  Vj  (0,  ac)  ,  djvj  (0,  k)  from  the  expansions 

N 

v(x,y,0,K)  =  Yxj  (0 ,K)(j)j(x,y), 

3=1 


N 

d]v  (x,  y,  0,/c)  =  X  (°>  K)  0j(x’  V%  7  =  1?  2. 

3=1 

Since  the  functions  v  (x,  y,  0,  ac)  ,  vz  (x,  y,  0,  /c)  are  known,  functions  Vj  (0,  /c)  ,  dzVj  (0,  ac)  can 
be  obtained  by  solving  a  system  of  linear  equations.  For  the  second  derivatives  d2Vj  (0,  k)  , 
we  assume  the  coefficient  c(x)  is  known  at  the  surface  of  measurements  2  =  0,  which  is  often 
the  case  in  practice.  So  the  second  derivatives  d 2Vj  (0,  /c)  can  be  computed  from  equation 
(2.22). 

To  approximate  numbers  in  (3.17),  we  combine  (3.15),  (3.16)  and  apply  the  least  squares 
minimization  in  L2,  i.e., 

min 


M7)  (z0)  C7J  {zq) 


(7) 


d]vj  (0,  ac) 


AC 


AC 


dt c,  7  =  0, 1, 2, 


(3.20) 
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where  (20)  :=  kj  (20) ,  t!p  (20)  :=  tj  (20).  By  imposing  the  critical  condition,  the  vector 
[k^  (20) ,  t(kf}  (^o)j,  the  minimizers  of  the  integral  (3.20)  is  the  solution  of  the  2x2  system 
of  linear  algebraic  equations 


where 


ak^  (z0)  +  (20)  = 

bkf]  (20)  +  ctj7)  (20)  = 


K 


kt 


d]Vj  (0,  k 


K 


K 


2  K 


tt2  ’ 


C  = 


—d,K, 

(3.21) 

-dn, 

(3.22) 

f  (Ik  1  1 

J  ft4  3k3  3F3 

The  determinant  of  the  system  (3.21),  (3.22)  does  not  vanish.  With  vectors  V  (20)  :  = 
(ki  (zo) , ...,  k^  (20))  /k,  V'  (20)  :=  (k[(z0) ,  ...,k'N  (z0)) /k  computed,  we  can  minimize  the 
functional  J\  (ai(/c))  for  the  hrst  layer.  Functions  p3  (21,  k),  p)  (21,  k)  and  p"  (zl7  k)  can 
be  evaluated  once  the  solution  is  obtained. 

Step  2.  Consider  the  second  layer  2  G  [21,22].  Denote  k^  (21)  =  kj  (21) ,  t^  (21)  = 
tj  (21) .  From  (3.15)  and  (3.16) 


d^Vj  (21,  k)  = 


k] 


(7) 


(*.)  ,  C(2.) 


K 


+ 


z  e  [Zi-i,  Zi) ,  K  E  [/c,  k],  7  =  0, 1,  2. 


The  same  procedure  can  be  applied  to  find  the  approximate  values  of  k ^  (21) ,  t'k’'1  (21),  i.e., 
we  minimize  the  least  squares  functional  in  L 2  (k,  k) , 


mm 


d'lvj  (21,  k)  - 


k) 


(7) 


(d) 


Ai) 


(d) 


~i  2 


K 


K~ 


(Ik 


and  obtain  a  system  similar  to  (3.12),  (3.22).  Once  the  numbers  k ^  (21)  ,t^  (21)  are  found, 
tails  V (21)  and  their  derivatives  V'{z\),  V"{z\)  are  calculated  from  (3.21),  (3.22)  as 


^  7  =  o,1,2. 

K 


Other  layers  are  treated  similarly. 
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4  Comparative  Analysis 

In  order  to  fully  investigate  the  numerical  performance  of  the  second  implementation  of  the 
convexification  algorithm,  we  conduct  in  this  section  a  comparative  study  on  ten  aspects 
of  this  algorithm  [4].  It  is  carried  out  with  three  configurations  for  the  coefficient  inverse 
problem.  Recall  that  we  consider  the  case  when  mine-like  inclusions  are  located  in  a  dry  soil 
(subsection  3.1).  The  simplest  case  is  with  the  homogeneous  background  medium  without 
mines.  One  realistic  configuration  is  with  a  single  mine,  the  one-mine  case.  The  other  is 
with  two  inclusions  and  we  consider  a  stone  and  a  mine  as  the  inclusions,  the  stone-mine 
case.  So,  our  result  below  show  that  we  can  differentiate  between  stones  and  mines  on  our 
images.  The  stone  and  mine  are  of  the  circular  shape  with  different  radius,  rstone  —  4.5cm 
and  rmine  =  5cm.  Assuming  that  the  stone  is  “more  wet”  than  the  dry  soil,  we  model  the 
stone  with  the  same  relative  dielectric  permittivity  as  the  wet  soil  with  5%  moisture  (Figure 
3),  thus,  esrtone  =  4.  For  the  case  with  a  single  mine,  the  center  of  the  mine  is  located  at 
the  point  Pm(3Qcrn,  7.5 cm),  whereas  for  the  stone-mine  case,  the  centers  of  the  mine  and  the 
stone  are  at  the  point  Pm (40cm,  7.5cm)  and  the  point  Ps(— 40cm,  7.5cm),  respectively. 

4.1  Efficiency 

All  our  computation  was  performed  on  a  workstation  with  the  CPU  of  AMD  Athlon(tm)  64 
Processor  3500+,  2.2  GHz  clock  speed,  4  GB  of  RAM,  512  KB  cache  size,  running  Novell 
Client  for  Linux  vl.O.  The  code  is  written  in  C,  optimized  up  to  three  levels,  and  the 
compiler  is  gcc-3.4.6.  To  demonstrate  the  efficiency  of  the  convexification  algorithm  with 
our  new  implementation,  we  consider  the  two  realistic  configurations.  For  the  case  with 
a  single  mine,  we  use  35  cubic  R-splines  for  functions  O-(x)  with  layer  size  Sz  =  0.005, 
k  G  [9, 10]  and  apply  Legendre  polynomials  of  degree  5  for  pseudo- frequency  approximation, 
and  add  2%  noise  in  the  input  data.  The  code  is  run  up  to  16  layers  and  the  iteration  stops 
when  the  gradient  of  the  objective  function  J  drops  below  the  threshold  5  =  0.01.  The 
computational  time  for  the  run  is  1  minute  and  14  seconds.  For  the  configuration  with  a 
stone  and  a  mine,  we  apply  39  cubic  R-splines  for  functions  (b3{x)  and  all  other  parameters 
share  their  corresponding  values  in  the  previous  case  with  a  single  mine.  It  takes  1  minute 
and  38  seconds  for  the  code  to  run  up  to  16  layers.  Thus,  the  total  depth  underneath  the 
ground  is  8  cm,  which  is  a  reasonable  depth  for  antipersonnel  land  mines. 

On  each  generic  Layer  $4,  the  error  of  the  minimizer  a  j  from  the  minimization  procedure 
is  bigger  than  that  with  a  tight  threshold,  e.g.,  5  =  0.001,  or  §  =  0.0001.  However,  the 
reconstructed  unknown  coefficient  er  for  each  configuration  is  still  good,  which  is  shown 
in  Figure  4.  We  should  point  out  that  although  our  computations  were  performed  on  a 
computer  with  a  large  memory,  i.e.,  4GB  of  RAM,  the  compiled  code,  which  is  rather  small, 
can  be  run  on  a  computer  with  normal  amount  of  memory,  e.g.,  512  Mb  of  RAM,  or  even 
smaller  size,  e.g.,  256  Mb  of  RAM.  The  memory  requirement  for  the  code  to  run  is  rather 
low.  The  numbers  of  degrees  of  freedom  for  the  two  cases  are  210  and  234,  respectively.  This 
points  towards  a  possibility  for  the  convexification  to  work  in  real  time  to  detect  and  image 
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Figure  4:  Reconstructed  er  for  the  efficiency  study  at  z  =4.75  cm,  top:  the  configuration 
with  a  single  mine,  bottom:  the  configuration  with  a  stone  and  a  mine. 


antipersonnel  land  mines  in  the  field. 

The  solution  of  the  forward  problem  has  a  singularity  at  the  point  Xq  =  (0,  — 10cm) , 
where  the  source  is  located.  This  influences  the  sensitivity  region  for  the  inverse  problem. 
Numerical  experiences  show  that  the  reconstructed  unknown  coefficient  er  is  not  sensitive 
near  the  singularity,  i.e. ,  near  the  center  of  the  domain.  We  have  also  investigated  the 
case  when  the  central  part  is  imaged  [5].  In  this  case  the  resulting  image  has  a  small 
dent  in  the  center  by  the  presence  of  the  singularity  near  the  source  location.  We  have 
explained  this  dent  in  [5].  Figure  5  represents  a  typical  image  when  the  central  part  is  in. 
Unless  otherwise  mentioned  explicitly,  each  graph  below,  is  the  cross-sectional  view  of  the 
reconstructed  unknown  coefficient  er  along  a:-axis  without  the  central  part  at  a  specific  value 
in  z  direction. 


4.2  Stopping  Criterion 

On  each  generic  Layer  #i,  the  steepest  descent  method  is  applied  to  find  the  unique  minimizer 
of  the  weighted  least  squares  objective  function  ,JX.  The  iterative  solver  has  the  form 


(n+i)  _  (n) 

a  i  ~  ai 


-  U’VJj 


An) 


where  a ^  >  0  is  the  step  size  at  ri-tli  iteration.  The  iteration  is  terminated  when  the 


gradient  of  the  function  drops  below  a  prescribed  threshold,  i.e., 


An) 


<  5.  We 
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Figure  5:  Reconstructed  coefficient  er,  the  case  with  a  stone  and  a  mine,  £=4.75  cm.  The 
dent  in  the  center  is  clearly  seen  and  was  explained  in  [5]. 

consider  three  thresholds  5  =  {1  x  1(U2,1  x  10-  "3, 1  x  10  4}  with  the  homogeneous  case 
and  run  the  code  np  to  16  layers.  The  lower  and  upper  limit  of  the  pseudo-frequency  is 
k0  =  9  and  R  =  10,  respectively,  i.e. ,  k  G  [9,10].  For  this  range  of  pseudo-frequency, 
Legendre  polynomials  of  degree  5  give  the  best  approximation,  which  is  used  as  the  basis 
for  ft  approximation  in  this  run.  Value  for  the  parameter  A  associated  with  the  Carleman 
Weight  Function  is  A  =  200.  We  use  39  splines  on  the  interval  x  G  [—A,  A]  =  [—0.6,  0.6] 
with  uniform  layer  size  8z  =  0.005.  Unless  otherwise  noted  explicitly,  this  uniform  layer  size 
has  been  used  throughout  the  subsequent  comparative  studies.  Figure  6  shows  the  result  of 
the  reconstruction  for  the  case  er  =  2.9,  i.e.,  in  the  absence  of  the  inclusion.  From  the  graph 
with  normal  scale,  the  relative  dielectric  permittivity  er  is  reconstructed  successfully  for  all 
stopping  criteria  though  with  very  small  deviations  from  the  exact  value.  The  difference  with 
the  three  stopping  criteria  is  self-evident  from  the  blow-up  view.  The  stringent  threshold 
5  =  1  x  10~4  does  not  give  us  a  better  accuracy.  Rather,  it  introduces  more  oscillations 
compared  with  the  relatively  loose  threshold  8  =  1  x  10~3. 

This  is  due  to  the  ill-posed  feature  of  coefficient  inverse  problems.  Accumulated  round-off 
error  also  contributes  this  since  more  iterations  are  involved  to  reach  the  stringent  criterion. 
It  is  a  delicate  issue  to  choose  an  appropriate  stopping  criterion.  If  it  is  too  loose,  i.e.,  5 
is  relatively  big,  then  the  minimizer  a*  from  the  steepest  descent  method  will  be  rather  far 
from  its  true  solution  a*.  This  will  result  in  degraded  accuracy  of  the  reconstructed  unknown 
coefficient,  which  is  shown  in  the  blow-up  view  for  the  threshold  8  =  1  x  10”2.  One  has  to 
take  into  account  the  error  with  spatial  approximation,  the  error  with  pseudo-frequency  ft 
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Figure  6:  Comparison  of  three  different  thresholds  at  z=8  cm,  top:  blow-up  view,  bottom: 
normal  scale  view. 

integration  and  layer  size  Sz  in  selection  of  a  proper  stopping  criterion. 

4.3  Layer  Size 

We  consider  three  different  layer  sizes  Sz  =  {2.5 x  10“3,  5xlCT3,  lx  10~2}  for  the  configuration 
with  no  inclusion,  i.e. ,  er  =  2.9.  The  other  parameters  are  the  same  as  in  the  previous  section 
on  stopping  criterion  with  8  =  1  x  10~3.  The  codes  run  up  to  32,  16  and  8  layers  for  the  three 
layer  sizes  considered,  respectively.  We  show  the  cross-sectional  view  of  the  reconstructed  er 
at  z  =  0.08  along  x-axis,  i.e.,  8  cm  underneath  the  ground  and  at  the  end  of  the  last  layer 
for  each  run.  The  comparison  is  shown  in  Figure  7.  On  the  graph  with  normal  scale,  the 
three  curves  are  close  to  each  other  and  er  is  reconstructed  very  accurately  for  each  layer 
size  considered.  From  the  blow-up  view,  one  can  appreciate  the  difference  -  the  thinner  the 
layer,  the  less  the  oscillations,  thus,  the  better  accuracy.  This  is  understandable  since  the 
quadratic  approximation  of  function  pj(z ,  hi)  with  respect  to  z  has  more  accuracy  with  layer 
size  decreased.  The  initial  values  of  a3+1,  b°+1,  c3+1  for  the  next  layer  also  have  less  errors 
with  smaller  Sz. 

4.4  Parameter  A 

We  study  the  impact  of  the  parameter  A  associated  with  the  Carleman  Weight  Function  on 
the  accuracy  of  the  reconstructed  er.  We  consider  A  =  {100,200,1000}  and  the  stopping 
criterion  8  —  1  x  10-3.  The  rest  of  the  parameters  share  the  same  values  as  those  in  the  section 
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Figure  7:  Comparison  of  three  different  layer  sizes  Sz  at  z  =8  cm,  top:  blow-up  view,  bottom: 
normal-scale  view. 

on  stopping  criterion.  Again,  we  work  with  the  homogeneous  situation  where  er  =  2.9.  The 
result  is  shown  in  Figure  8.  The  relative  dielectric  permittivity  er  is  reconstructed  with 
a  good  accuracy  for  all  values  of  A  considered.  The  difference  shows  up  on  the  blow-up 
view.  Smaller  value  of  A  generates  more  oscillations  and  larger  value  of  A  tends  to  create 
less  oscillations.  This  does  not  necessarily  mean  that  the  larger  value  of  A  will  have  a  better 
accuracy  for  the  reconstructed  er.  The  overall  performance  for  A  =  200  seems  to  be  better 
than  that  for  A  =  1000.  The  goal  of  the  introduction  of  the  Carleman  Weight  Function 
is  to  ensure  the  weighted  least  squares  objective  function  J\  to  be  strictly  convex  on  each 
generic  layer.  Certainly  the  parameter  A  has  impact  on  the  solution  a*  from  the  minimization 
procedure,  thus  affecting  the  reconstructed  er. 

From  the  results  shown  in  Figure  8,  we  conclude  that  the  objective  function  Jlx  has  a 
unique  minimizer  when  parameter  A  varies  in  a  large  range.  Furthermore,  the  results  conform 
to  the  error  estimate  (2.34),  which  means  that  the  algorithm  is  very  robust.  With  all  other 
parameters  fixed,  there  is  an  unique,  optimal  value  of  the  parameter  A*.  Actually,  we  have 
employed  the  implicit  rule  A  x  Sz  =  1,  i.e. ,  taking  the  inverse  of  layer  size  Sz  as  its  value 
for  this  parameter  in  the  rest  of  our  comparative  study.  Computational  results  for  various 
realistic  configurations  show  that  this  choice  works  quite  well. 

4.5  Dimension  of  the  Basis  for  k  Approximation 

An  orthonormal  basis  on  the  interval  [ko,  k]  has  been  applied  to  approximate  functions 
p (z,k)  and  p,(^,K)  with  respect  to  pseudo- frequency  k.  We  choose  Legendre  polynomials 
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Figure  8:  Comparison  of  three  different  values  for  the  parameter  A  at  £=8  cm,  top:  blow-up 
view,  bottom:  normal-scale  view. 


up  to  degree  K  as  the  basis  {/*  (k)}^0.  One  advantage  of  introducing  such  a  basis  is  that  all 
the  terms  involving  ^-integration  in  the  weighted  least  squares  objective  function  Jlx  can  be 
computed  explicitly  in  contrast  to  the  composite  Simpson’s  or  trapezoidal  rules.  Functions 
of  the  orthonormal  basis  are  defined  as 


li^K)  . 


2i  +  1  f  2k  —  (k0  +  k) 


K  —  Kq 


K 


K0 


where  Lj(£)  is  the  classical  Legendre  polynomial  of  degree  i  defined  on  [-1,1]  and  has  the 
property 
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The  orthonormal  basis  share  the  property 


K, 

J  k(K)lj(K)dK  =  Sij , 

where  8^  is  the  Kronecker  delta.  On  the  pseudo- frequency  interval  k  G  [k0,k]  =  [9,10], 
Legendre  polynomials  of  degree  5  have  the  best  approximation  of  functions  p(0,  k)  and 
p2(0,fi:)  compared  with  other  degrees.  Indeed  for  degree  5  polynomial  approximation, 
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the  maximum  absolute  and  maximum  relative  errors  are  very  small  for  both  functions: 
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ire/  oo (0;  /c)  =  3.62  x  10~4.  We  want  to  investigate  if  there  is  any  pronounced  difference 
in  the  reconstructed  unknown  er  with  these  two  polynomial  approximations.  For  the  com¬ 
parison,  we  consider  one  realistic  configuration,  the  one-mine  case.  The  center  of  the  mine 
is  located  at  the  point  Pm(0. 3m,  0.075m)  =  Pm (30cm,  7.5cm).  We  use  43  splines  for  spatial 
approximation  and  the  threshold  6  =  0.001  for  the  stopping  criterion. 

Figure  9  shows  the  comparison  result.  The  mine  and  the  background  medium  have  been 
clearly  identified  for  both  polynomial  approximations.  On  the  graph  with  normal  scale,  the 
two  curves  are  very  close  to  each  other.  From  the  two  blow-up  views,  one  can  recognize  the 
difference  though  not  too  much.  Degree  5  gives  relatively  better  overall  performance  than 
degree  2.  Compared  with  the  big,  several  orders’  difference  in  approximation  of  functions 
p(0,  k)  and  p,(0,  k),  the  discrepancy  of  er  is  rather  small  with  the  two  different  polynomial 
degrees.  This  is  due  to  the  fact  that  the  tail  function  x(x,  z )  dominates  the  contribution 
for  the  function  v(x,z,k),  and  therefore  we  could  take  a  very  short  interval  k  G  [9,10]  for 
the  pseudo-frequency  integration.  On  this  interval  the  approximation  errors  with  Legendre 
polynomials  of  degree  2  and  5  are  all  very  small.  Thus,  this  small  error  propagates  through 
the  minimization  procedure  and  exerts  a  slight  influence  on  the  final  reconstructed  unknown 
coefficient. 


4.6  Noisy  Data 

We  consider  the  perturbed  input  data  by  adding  different  noise  levels  to  the  input  data 
(p2(x,  K)  and  "02 (^5  K).  I n  principle,  we  should  add  noise  to  the  original  data  <p  and  ^  in  (3.6) 
and  then  use  a  regularization  procedure  to  differentiate  functions  tp1  and  rj)1  in  (2.14)  with 
respect  to  k.  Regularization  can  be  done,  for  example  similarly  to  subsection  6.1  of  [34], 
However,  to  avoid  additional  complications,  instead  of  considering  (p(x,  k)  and  ^(x,  k),  in  this 
study  we  introduce  multiplicative  noise  in  the  input  data  tp2(x,hi)  and  i/j2(x,k).  For  every 
discrete  value  aq  G  [—4,  A]  and  Kj  G  [tt0,  k],  each  of  the  function  (p2 (xi:  Kj)  and  Kj)  is 

a  matrix.  The  elements  in  the  two  matrices  are  perturbed  by 

Xp2  (XU  Kj)  =  (p2  (Xi,  Kj)  (1  +  Cij)  ,  -02  (xi,  fy)  =  i>2  (xi,  Kj)  (l  +  T}^)  , 

where  and  T)tJ  are  normally  distributed  pseudo-random  numbers  with  zero  means.  The 
standard  deviations  are  selected  to  create  noise  levels  of  2%,  5%  and  10%  for  the  input  data 
<P2(x,k)  and  'i/j2(x,k).  The  sample  number  is  taken  to  be  200  for  each  noise  level.  The 
parameters  for  this  setup  are  the  same  as  those  in  the  previous  subsection  with  polynomials 
of  degree  5  for  pseudo-frequency  approximation.  The  results  with  the  perturbed,  noisy  input 
data  are  shown  in  Figure  10. 

For  the  reconstructed  unknown  coefficient  sr,  higher  noise  level  creates  more  oscillations 
and  there  is  small  undershoots  near  the  mine  with  the  noise  level  of  10%.  Nevertheless, 
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Figure  9:  Comparison  of  polynomial  approximations  with  degree  2  and  5  at  z— 5.5  cm,  top: 
normal-scale  view,  middle:  blow-up  view  of  the  upper  portion,  bottom:  blow  up  view  near 
the  mine. 
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Figure  10:  Comparison  of  three  different  noise  levels  at  z=5.5  cm,  top:  normal-scale  view, 
bottom:  blow-up  view  near  the  mine. 
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Figure  11:  Comparison  between  the  case  with  the  tail  and  the  tail-free  case  at  z=4.75  cm, 
top:  with  tail,  bottom:  tail-free. 


the  overall  performance  with  all  the  noise  levels  are  still  good  and,  both  the  mine  and  the 
background  medium  have  been  clearly  identified  in  terms  of  their  material  property  er. 

4.7  Tail  V(x) 

As  we  have  pointed  out  before,  the  tail  function  V(x)  dominates  the  contribution  in  the 
pseudo-frequency  k  integration.  We  compare  the  reconstructed  unknown  coefficient  er  for 
the  case  with  tails  and  the  case  without  tails,  the  tail-free  case  and  consider  the  stone-mine 
realistic  configuration.  The  center  of  the  stone  is  located  at  the  point  Ps{~ 0.4m,  0.075m)  = 
Ps{— 40cm,  7.5cm),  i.e. ,  the  same  depth  7.5  cm  underneath  the  ground  as  the  mine,  whose 
center  is  at  the  point  Pm(0. 4m,  0.075m)  =  Pm (40cm,  7.5cm).  The  ^-interval  for  the  inverse 
problem  is  :=  [—A,  A]  =  [—0.7,  0.7].  Each  computation  is  performed  with  51  splines  with 
uniform  layer  size  8z  =  0.005m  =  0.5cm  and  6  =  1  x  10”4  as  the  threshold  for  the  stopping 
criterion.  The  lower  limit  for  k  integration  is  Ko  =  8  and  degree  6  for  Legendre  polynomial 
basis.  The  comparison  is  presented  in  Figure  11.  The  tail-free  case  essentially  gives  us 
nothing  on  the  reconstructed  er,  i.e.,  it  could  not  determine  the  material  property  (er)  of  the 
heterogeneous  media,  let  alone  distinguish  the  stone  and  the  mine  from  their  background 
medium.  Yet,  the  variation  of  er  for  the  tail-free  case  still  indicates  the  locations  of  the  stone 
and  the  mine. 
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Figure  12:  Comparison  of  four  different  number  of  splines  at  z=5.5  cm,  top:  normal-scale 
view,  second  row:  blow  up  view  near  the  mine,  third  row:  blow-np  view  of  the  left  portion, 
bottom:  blow-np  view  for  x  G  [0.32,  0.44]  . 

4.8  Dimension  of  the  Basis  for  x  Approximation 

To  approximate  function  q(x,  z,  k)  and  its  partial  derivatives  with  respect  to  spatial  variable 
x,  we  introduce  cubic  5-splines:  (x)  ,  ■  Each  spline  function  is  piecewise  defined 

on  a  subset  of  the  entire  interval  [a,  b],  yet  globally  the  interpolating  function  /(x)  has  the 
property  f(x)  €  C2[a,b\.  It  is  well  known  that  it  has  the  desirable  approximation  property, 
i.e. ,  if  /  G  C4[a,b],  then 

11/  -/Hoc  <  ^||/(4)|U4, 

where  h  is  the  distance  between  two  consecutive  evenly-spaced  knots.  We  study  the  effect 
of  the  dimension  N  of  the  basis  on  the  resolution  of  the  final  reconstructed  er,  i.e.,  whether 
using  more  splines  will  give  us  better  accuracy  of  £r.  The  parameters  share  the  same  values 
as  in  Subsection  4.5  and  for  k  approximation,  we  choose  K  =  5  and  k  G  [9, 10].  The  results 
are  shown  on  Figure  12.  On  the  graph  with  the  normal  scale,  all  the  4  curves  are  close  to 
each  other,  but  from  the  3  blow-up  views,  one  can  identify  the  difference.  Near  the  mine, 
there  is  small  undershoot  of  £r  for  N  =  {43,47}.  In  this  local  region,  the  best  performance 
seems  to  be  the  case  with  N  =  39. 

At  the  left  portion  of  the  graph,  there  are  mild  oscillations  for  N  =  {43,47},  and  the 
oscillation  for  N  =  35  is  rather  small.  The  performance  with  N  =  39  is  in-between  the 
performances  for  N  =  35  and  N  =  43.  The  best  resolution  of  er  in  this  local  area  is  for 
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the  case  with  N  =  35.  From  the  blow-up  view  for  x  G  [0.32,  0.44]  and  near  the  interface  of 
singularity,  overshoots  exist  for  all  the  cases  with  the  four  different  splines  numbers.  The 
curve  for  N  =  35  shows  the  maximum  overshoot  while  the  minimum  overshoot  exhibits 
on  the  curve  for  N  =  47.  The  best  overall  performance  in  this  short  interval  seems  to 
be  for  N  =  39.  Thus,  on  the  whole  interval  x  G  [—0.6,  0.6],  the  highest  accuracy  of  the 
reconstructed  er  is  for  the  case  with  N  =  39.  The  main  reason  why  further  more  splines 
does  not  improve  the  accuracy  of  the  reconstructed  er  is  due  to  the  inherent  ill-posed  nature 
of  coefficient  inverse  problems.  Accumulated  round-off  error  also  contributes  but  to  a  less 
degree.  Another  disadvantage  with  more  number  of  spline  approximation  is  the  increased 
computational  load,  i.e.,  more  running  time  and  more  required  memory  of  the  machine.  This 
might  be  a  problem  for  detecting  the  land  mines  in  real  time. 


4.9  Lower  Limit  kq 


Once  R  is  decided  using  the  asymptotic  behavior  of  function  w(x,  s)  as  s  — >  oo,  there  is  an 
arbitrariness  in  selecting  the  lower  limit  k,q  for  pseudo-frequency  integration.  Theoretically, 
fto  should  be  as  small  as  possible  in  order  to  have  better  accuracy  for  k  integration.  The 
problem  with  smaller  value  of  kq  is  that  the  dimension  of  the  Legendre  polynomial  basis  will 
also  have  to  be  correspondingly  higher.  For  example,  if  we  take  Kq  =  1,  Legendre  polynomials 
of  degree  17  give  the  best  approximation  of  functions  p(0,  k)  and  pz(0,  k).  This  is  in  a  sharp 
contrast  with  the  case  for  Ko  =  8,  where  the  degree  6  gives  the  best  approximation  for 
both  functions.  Here  we  want  to  study  the  impact  of  different  values  of  /c0,  i.e.,  k,o  =  8 
and  Ko  =  9  on  the  resolution  of  our  reconstructed  £r.  We  record  the  maximum  absolute 
and  maximum  relative  errors  with  degree  6  polynomial  approximations  for  both  functions: 
MEpO.O  =  l-«  x  10-®,  lipISfO.K)  =  3.47  x  KT7;  felELM  =  2.45  x  KT9, 
I'Sp.lSto  (0,  k)  =  3.29  x  10  7 .  Each  error  is  on  the  same  order  as  its  corresponding  one  with 
degree  5  approximation  for  k0  =  9,  though  a  little  larger.  For  this  study,  we  consider  a 
realistic  configuration,  the  stone-mine  case.  The  rest  of  parameters  are  the  same  as  those  in 
subsection  4.7.  Figure  13  shows  the  results  for  the  comparison. 

The  two  curves  are  almost  inappreciable  from  the  graph  with  normal  scale,  and  there  is 
no  overshoot  near  the  stone,  nor  undershoot  near  the  mine.  The  stone,  the  mine  and  the 
background  medium  have  been  correctly  and  sharply  identified.  The  difference  shows  up  on 
both  blow-up  views.  The  performance  of  the  reconstructed  er  is  a  little  better  with  kq  =  8. 
This  is  due  to  the  extra  information  with  k  integration  for  k  G  [8,9],  in  the  presence  with 
the  case  for  k0  =  8  and  lost  in  the  case  for  k0  =  9. 


4.10  Initial  Guess  a? 

To  start  the  layer  stripping  procedure,  we  need  to  have  the  initial  guess  a]1  for  the  first 
layer  z  G  [0,5z].  There  are  two  ways  to  obtain  a,.  The  first  approach  needs  one  more 
datum  wzz(x,  0,  k)  besides  the  two  given  Cauchy  data,  w(x,  0,  k)  and  wz(x,  0,  k).  This  extra 
datum  can  be  obtained  by  solving  the  forward  problem.  In  order  to  distinguish  from  the 


36 


Figure  13:  Comparison  of  two  different  values  of  Koat  £=4.75  cm,  top:  normal-scale  view, 
middle:  blow  up  view  near  the  stone,  bottom:  blow-up  view  near  the  mine. 
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other  method,  we  term  this  approach  as  “External  data”  for  the  obvious  reason.  The  second 
method  is  by  utilizing  equation  (2.27).  For  the  first  layer  and  at  z  —  0,  by  the  quadratic 
approximation  of  p (z,k)  in  (2.27),  we  have  p(0,  k)  =  cx,  p'(0,  k)  =  bx,  and  p"(0,  k)  =  ax. 
Plugging  these  expressions  into  equation  (2.27)  and  re-arranging  the  second  term  on  the  left 
hand  side  yields 


B ax=f  Ci,  bx,  /  c1(z0,T)dr,  /  bx  (z0,  r)  dr,  x(z0),  x'(*o),  k 


The  vector-valued  function  f  can  be  evaluated  since  all  its  arguments  are  known.  As  matrix 
B  is  non-singular,  we  could  solve  the  above  linear  equation  for  ax  to  get  the  initial  guess  a)1 
for  the  first  layer.  We  name  this  method  “By  formula”.  For  this  run,  the  parameters  take 
the  same  values  as  their  corresponding  ones  in  subsection  4.7  for  the  case  with  tail.  The 
comparison  results  are  shown  in  Figure  14.  On  the  graph  with  normal  scale,  the  two  curves 
are  nearly  indistinguishable. 

From  the  blow-up  view,  one  can  tell  the  difference:  the  by-formula  method  is  a  little  better 
than  the  approach  with  external  data  in  terms  of  resolution  of  the  reconstructed  unknown 
coefficient  er.  The  convergence  history  for  Layer  #10  shows  that  the  by-formula  method 
converges  a  little  slower  than  the  approach  with  external  data.  The  former  needs  33  iterations 
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compared  with  29  iterations  for  the  latter  to  reach  the  stopping  criterion 
The  iteration  number  between  the  two  approaches  differs  only  slightly.  We  consider  this 
comparison  is  a  testimony  of  the  fact  that  the  convexihcation  algorithm  is  indeed  a  globally 
convergent  method. 


5  Summary  and  Conclusions  of  Comparative  Studies 

A  systematic  comparative  study  of  the  globally  convergent  convexihcation  algorithm  has 
been  carried  out.  Below  is  a  summary  of  our  findings. 

•  An  appropriate  stopping  criterion  for  the  iteration  of  gradient-like  method  is  central 
to  the  accuracy  of  the  reconstructed  unknown  coefficient.  A  stringent  criterion  does 
not  necessarily  guarantee  better  accuracy.  To  choose  a  proper  threshold,  one  needs 
to  consider  the  errors  of  spatial  approximation,  of  pseudo-frequency  integration  and 
layer  size.  With  a  known  configuration  and  by  trial  and  error,  one  can  calibrate  the 
threshold  and  apply  it  to  other  unknown  realistic  situations.  Our  numerical  studies 
show  that  for  realistic  configurations,  the  difference  of  the  reconstructed  unknown 
coefficient  er  between  a  stringent  and  a  loose  threshold  is  small.  On  the  other  hand, 
the  stringent  threshold  requires  much  more  computational  time.  This  means  that  with 
a  loose  threshold,  we  can  still  generate  a  good  image  of  the  unknown  coefficient,  and 
more  importantly,  the  efficiency  study  shows  that  it  is  possible  to  run  the  algorithm  in 
real  time  for  identification  and  imaging  of  antipersonnel  land  mines  in  the  field. 
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Figure  14:  Comparison  of  two  different  initial  guesses  of  at  z=4.75  cm,  top:  normal-scale 
view,  middle:  blow-up  near  the  mine,  bottom:  convergence  history  for  layer  #10. 


39 


•  Thinner  layer  size  has  better  resolution  of  the  reconstructed  unknown  coefficient  er, 
which  provides  a  numerical  confirmation  of  the  estimate  (2.35). 

•  The  weighted  least  squares  objective  function  is  strictly  convex  for  the  parameter  A 
associated  with  the  Carleman  Weight  Function  varies  in  a  large  range.  As  long  as 
A  varies  in  this  range,  its  influence  on  the  accuracy  of  the  reconstructed  er  is  not 
pronounced,  which  confirms  the  error  estimate  (2.34).  The  inverse  of  the  layer  size  can 
be  a  good  choice  for  the  value  of  this  parameter. 

•  A  basis  of  low  dimension  can  be  applied  for  the  pseudo-frequency  approximation  and 
generate  good  accuracy  on  er.  This  is  due  to  the  dominance  of  the  tails  in  the  pseudo¬ 
frequency  integration. 

•  A  higher  noise  level  affects  more  severely  the  quality  of  the  reconstruction  of  the  un¬ 
known  coefficient  er.  Up  to  the  10%  noise  level,  the  quality  of  the  reconstructed  er  is 
still  good. 

•  For  each  specific  problem,  there  is  an  appropriate  dimension  of  the  basis  for  spatial 
x  approximation.  Due  to  the  inherent  ill-posed  nature  of  coefficient  inverse  problems, 
basis  with  even  higher  dimension  does  not  further  improve  quality  of  the  reconstructed 
unknown  coefficient. 

•  The  lower  bound  of  the  pseudo-frequency  may  be  selected  quite  close  to  the  upper 
bound  due  to  the  fact  that  tails  dominate  pseudo- frequency  integration.  Small  change 
of  the  lower  limit  has  negligible  impact  on  the  reconstructed  unknown  coefficient. 

•  Fitting  the  correct  tails  is  crucial  for  the  algorithm  to  work  properly  and  efficiently. 

•  The  initial  guesses  obtained  by  the  methods  of  “External  data”  and  “By  formula”  differ 
only  slightly  in  terms  of  resolution  of  reconstructed  unknown  coefficient  and  speed  of 
convergence.  This  conforms  well  with  the  theoretically  established  strict  convexity  of 
the  functional  J\. 

•  Overall,  the  comparative  study  demonstrates  the  robustness  of  the  convexification  algo¬ 
rithm.  Efficiency  with  the  new  implementation  shows  the  possibility  for  the  algorithm 
to  be  applied  in  real  time  to  detect  and  image  antipersonnel  land  mines  in  the  field. 

The  above  comparative  analysis  provides  valuable  references  to  further  development  of 
the  convexification  algorithm  for  a  broad  class  of  CIPs. 

6  Convexification  Method  for  an  Inverse  Problem  for 
an  Elliptic  Equation 

In  [7]  we  have  also  considered  the  inverse  problem  for  the  elliptic  equation  (6.1)  in  the  case 
when  the  running  pseudo  frequency  of  the  above  case  is  replaced  by  the  running  source 
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Xs.  The  source  is  running  along  a  line.  This  corresponds  to  the  constant  current  in  the 
case  of  search  of  e.g.,  land  mines  and  underground  bunkers  using  the  method  of  Electrical 
Impedance  Tomography  (E1T).  The  soil  is  probed  by  the  constant  current  then.  In  the  case 
of  EIT  the  function  u(x,xs)  in  (6.1)  is  the  voltage  at  the  point  x  generated  by  the  constant 
current  at  the  source  point  Xs.  In  the  case  of  optical  imaging  of  diffuse  media  this  corresponds 
to  the  so-called  Constant  Wave  (CW)  light.  The  optical  imaging  through  the  diffuse  media 
can  be  applied  to  image  targets  on  battlefields  through  smogs  and  flames,  as  well  as  to 
medical  optical  imaging.  Interestingly,  the  diffuse-like  propagation  of  light  in  the  battlefield 
application  would  even  be  helpful,  because  even  if  the  direct  laser  beam  would  “miss”  the 
target,  one  might  still  image  it  because  photons  would  still  “sense”  that  target  due  to  the 
diffusion  of  light.  A  simplified  version  of  the  convexification  for  this  inverse  problem  was 
considered  in  [9]  and  was  reported  in  the  annual  2006  report  on  this  project  [21]. 

The  reason  why  the  above  scheme  of  the  convexification  works  well  for  the  case  of  fre¬ 
quency/time  dependent  data  is  that  there  exists  a  “proper”  asymptotic  behavior  (2.11)  of 
the  solution  of  the  corresponding  elliptic  equation  when  the  pseudo-frequency  tends  to  the 
infinity  in  this  case.  However,  in  the  case  of  the  source  dependent  data  a  “proper”  asymp¬ 
totic  behavior  as  the  source  position  tends  to  the  infinity  is  unknown.  This  means  that  we 
need  a  special  treatment  of  the  so-called  “tail” .  The  tail  is  unknown  a  priori  and  appears 
due  to  the  truncation  of  an  improper  integral  with  respect  to  the  source  position,  see  (6.11). 
In  order  to  apply  our  layer  stripping  procedure,  we  need  to  approximate  the  tail.  Thus,  we 
use  a  heuristic  approach  of  approximating  the  tail.  At  the  same  time,  we  point  out  that  if 
the  tail  is  given,  then  the  global  convergence  of  our  layer  stripping  procedure  is  rigorously 
guaranteed. 

6.1  Statement  of  the  Inverse  Problem  and  Applications 

6.1.1  The  Inverse  Problem 

For  the  sake  of  generality  we  consider  the  3-D  case  with  x  =  (x\,X2,  z )  el3.  The  2-D  case, 
for  which  our  numerical  experiments  are  conducted,  is  both  simpler  and  similar.  Let  the 
function  a(x )  >  const.  >  0,a  6  C1  (M3)  and  a(x)  =  k2  =  const.  >  0  for  x  G  (|x|  >  R} , 
where  R  is  a  positive  number.  Let  the  function  u(x,s),  where  s  is  a  parameter,  satisfies  the 
following  elliptic  equation 

A xu  —  a(x)u  =  —5  (x  —  xs)  in  M3  (6.1) 

with  the  conventional  condition  at  the  infinity 

lim  u(x,  s)  =  0,  Vs  E  M.  (6.2) 

|#|— >oo 

The  source  position  is  at  the  point  Xs  and  when  s  is  changing,  the  source  is  running  along 
the  line  l  =  {x  =  (aq,  0,  zm)} ,  where  zm  =  const.  >  0.  We  also  assume  that 

a(x)  =  k2  in  le,  (6.3) 
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where  h  is  a  small  neighborhood  of  the  line  /. 

If  the  function  a(x)  =  A;2,  then  the  fundamental  solution  of  the  equation  (6.1)  with  the 
condition  (6.2)  is 


u0(x,s) 


exp  (—k  |  a;  —  xs|) 
47T  \x  —  xs\ 


Hence,  we  seek  the  solution  of  the  problem  (6.1),  (6.2)  in  the  form  u  =  Uq  +  u,  where  the 
function  u  E  C2  (M3)  and  satishes  the  following  conditions 


A xu  —  a(x)u  =  (a(x)  —  k2)  u0  in  M3, 


(6.4) 


lim  u(x,xs)  =  0.  (6.5) 

|  x  |  — >oo 

Uniqueness  and  existence  results  of  the  problem  (6.1),  (6.2)  for  functions 

u  E  C2+a  (M3\{|a;  —  x0|  <  e})  follow  from  the  classic  theory  of  elliptic  equations,  see, 
e.g.,  [].  Here  £  >  0  is  an  arbitrary  number,  a  E  (0, 1)  and  C2+a  are  Holder  spaces. 

Let  C  M3  be  a  rectangular  prism 


H  =  {x  :  -A  <  xi,x2  <  A,  z  E  (z0,  zm  -  6*)}  , 

where  A  and  L  are  positive  numbers.  Denote  T  =  n  {z  =  ^o}-  We  are  not  setting  9  :=  0 
because  we  want  to  avoid  working  with  the  line  /,  where  the  singularities  of  the  function  u 
occur.  We  consider  the  following 

Inverse  Problem.  Suppose  that  the  function  a(x)  is  unknown  inside  of  the  domain  Q, 
known  outside  of  it,  and  it  is  also  known  in  H  fl  {z  E  (zm  —  9,  zm )}  .  Determine  the  function 
a(x)  in  0n{z6  (zo,  zm  —  9)}  given  the  following  functions  cp  and  ^ 

u  (x,  s)  —  <p  (x,  s ) ,  uz  (x,  s)  =  ip  (x,  s ) ,  for  x  E  T,  s  E  [so,  s]  ,  (6.6) 

where  So  and  s  are  two  numbers.  Thus,  functions  <p  and  pj  represent  the  Dirichlet  and 
Neumann  data  caused  by  the  transmitted  signal. 


6.1.2  Applications 

In  applications  discussed  below  it  would  be  more  natural  to  consider  equation  (6.1)  either 
in  a  bounded  domain  or  in  the  half-space  {z  <  zm } .  In  principle,  these  cases  can  also  be 
incorporated  in  the  scheme  of  the  convexihcation.  However,  this  would  create  some  additional 
difficulties  in  the  forward  problem  solution  near  the  boundary.  Thus,  we  consider  the  forward 
problem  (6.1),  (6.2)  in  the  entire  space,  for  brevity. 

1.  Electrical  Impedance  Tomography 

One  of  applications  of  the  EIT  is  in  the  search  of  land  mines  and  underground  bunkers 
via  probing  the  ground  by  the  constant  current  at  different  source  locations.  Let  v(x,  s )  be 
the  voltage  generated  by  the  source  of  the  constant  current  located  at  (s,  0,  zm)  and  let  a(x) 
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be  the  electric  conductivity  of  the  medium,  a(x)  >  const.  >  0.  Then  the  function  v(x,  s ) 
satisfies  the  following  equation 

V  •  (a(x) Vv)  =  -6  (xi  -  s,  x2,  z  -  zm) 


and 


lim  v(x,  s )  =  0. 

|tc|  — >00 


Replacing  the  function  v  with  the  function  u  =  vffo  and  assuming  that  o  (xi,  0,  zm) 
obtain  equation  (6.1)  with 

a(x)  = 


1,  we 
(6.7) 


Hence,  assuming  that  the  function  <j(x)  is  known  near  the  surface  T,  we  arrive  at  the  above 
inverse  problem  with  the  unknown  coefficient  in  the  form  (6.7). 

2.  Optical  Diffusion  Tomography 

In  Optical  Diffusion  Tomography,  lasers  with  the  CW  (constant  wave)  light  are  used 
as  light  sources.  The  first  application  of  the  Optical  Diffusion  Tomography  is  in  optical 
medical  imaging  of  tumor-like  abnormalities  both  in  human  organs  and  small  animals  using 
Near  Infrared  (NIR)  light  with  the  wavelength  of  light  somewhere  between  500  and  1000 
nanometers.  The  second  feasible  application  is  in  optical  imaging  of  targets  on  battlefields 
via  smogs  and  flames.  Both  cases  of  transmitted  and  back  reflected  light  are  feasible  for 
both  applications.  The  light  source  should  move  along  a  straight  line  and  the  output  light 
should  be  measured  at  a  part  of  a  surface. 

Let  u(x,xs)  be  the  light  intensity  at  the  point  x  due  to  the  light  source  located  at  the 
point  xs.  It  can  be  derived  from  the  well  known  literature  sources  (see,  e.g.,  [31])  that  the 
function  u(x,xs)  is  the  solution  of  the  problem  (6.1),  (6.2)  with  the  coefficient  a  (x)  as 


a(x)  =  3  (/!>„)  (x) , 


(6.8) 


where  ffs  is  the  reduced  scattering  coefficient  and  fia  is  the  absorption  coefficient  of  the 
medium.  Both  coefficients  are  measured  in  [1/cm].  The  reduced  scattering  coefficient  p!s  is 
assumed  constant  here.  This  is  reasonable  because  in  NIR  applications  the  coefficient  ffs 
usually  changes  quite  slowly  with  respect  to  x  for  this  spectrum  of  light  waves,  whereas  the 
coefficient  )ia  changes  significantly.  Furthermore,  pa  can  be  used  for  the  diagnostics.  In  the 
case  of  imaging  of  targets  through  flames  these  targets  are  usually  impenetrable  for  light, 
meaning  that  fia  —  oo  in  them.  However,  one  can  model  those  targets  as  ones  with  finite, 
though  large  values  of  fia,  i.e.,  the  ones  having  large  contrasts  with  the  surrounding. 


6.2  The  Convexification 

Let  xs  =  (s,  0,  zm) ,  where  the  running  parameter  s  characterizes  the  changing  source  position 
of  the  source.  Below  we  consider  the  function  u(x,  s)  and  all  related  functions  only  in  the 
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above  domain  of  interest  fl  By  the  maximum  principle  u(x,  s )  >  0.  Hence,  similarly  with 
section  2  we  consider  the  function 


v(x,  s )  =  1iim(x,  s). 


By  (6.1) 


Av  +  (Vv)2  =  a(x)  in  fh 


(6.9) 


To  remove  the  unknown  coefficient  a(x)  from  (6.9),  differentiate  (6.9)  with  respect  to  s  and 
denote  g(x,s )  =  dsv(x,s).  Hence, 


A q  +  2 Vq  •  Vv  =  0. 


(6.10) 


We  now  need  to  express  the  function  v  via  the  function  q.  This  expression  is  obviously  given 
by 

S 

«(*,.)  =  -/9(*,T)rfr +  »(*,*)■  (Ml) 

s 

We  set  s  to  be  a  large  number.  We  call  the  function  v  (x,  s)  “tail” .  In  the  frequency  dependent 
case  the  function  v  (. x ,  s)  was  dropped,  for  s  »  1  because  v  (x,  s)  ~  0  in  that  case.  However, 
the  latter  is  not  true  in  our  case.  Hence,  we  need  a  special  treatment  to  approximate  the 
tail.  We  focus  now  on  an  approximation  of  the  function  q.  Substituting  (6.11)  in  (6.9),  we 
obtain 

S 

Aq-2Vq-  j  q(x,  t)cIt  +  2Vq  ■  Vv  (x,  s)  —  0.  (6.12) 

S 

Also,  conditions  (6.6)  imply  that 

q(x,s)  =  tp1(x,s),qz  (x,s)  =  (x,  s) ,  for  x  e  T,  s  e  [s0,  s\ ,  (6.13) 


where 


ib 

Vi{x,s)  =  dshup(x,s),'iljl  (x,s)  =  —(x,s)  - 


~2^ 


(x,s) 


Although  the  calculation  of  the  derivative  with  respect  to  s  is  an  ill-posed  procedure,  but  it 
can  be  handled  via  a  regularization  method,  see,  e.g.,  [34]  for  a  simple  method. 

We  have  obtained  the  Cauchy  problem  for  the  nonlinear  integral-differential  elliptic  Par¬ 
tial  Differential  Equation  (6.12)  with  the  Cauchy  data  (6.13).  Its  numerical  treatment  was 
similar  with  the  one  described  in  section  2.  As  to  the  basis  functions  in  (2. 21), we  have  chosen 
trigonometric  functions.  As  to  the  tail,  we  have  used  the  asymptotic  formula  for  the  solution 
of  the  forward  problem  (6.1),  (6.2)  and  the  available  data  to  get  an  approximate  value  of 
the  tail.  Interestingly,  we  have  discovered  that  it  is  sufficient  to  have  only  three  (3)  source 
positions  for  the  tail  and  only  two  (2)  source  position  for  the  subsequent  reconstruction  via 
the  convexihcation. 


44 


6.3  Some  details  of  numerical  experiments 

We  now  describe  some  ideas  of  numerical  experiments.  For  the  forward  problem,  we  calculate 
the  solution  of  the  diffusion  equation 

A u  —  a(x,  z)u  =  —5  (x  —  s,  z  —  zm)  (6.14) 

with  the  conventional  condition  at  the  infinity 

lim  u(x,z,s)=  0,  (6.15) 

|(x,z)|-KX> 

where  the  physics  of  the  function  a(x,z )  is  defined  in  (6.8).  We  have  considered  a  medical 
application.  However,  the  above  application  to  imaging  of  targets  on  battlefields  can  also  be 
considered  after  a  proper  re-scaling.  Consider  the  rectangle  12, 

12  =  {(x,  z)  :  5 cm  <  x  <  15 cm,  5cm  <  z  <  10cm}  . 


We  assume  that 

a(x,  z)  =  k2  =  const.  >  0  in  M2\12.  (6.16) 

We  assume  that  in  (6.14)  the  source  position  (s,  zrn)  is  running  along  the  right  side  of  0,  i.e., 
zm  =  L  =  10cm.  Also,  consider  a  bigger  rectangle 

Qo  =  {(a;,  z)  :  0 cm  <  x  <  20 cm,  0 cm  <  z  <  15 cm}. 

The  reason  why  we  consider  the  rectangle  120  along  with  the  rectangle  12  is  that  it  is  natural 
to  approximate  the  solution  of  the  problem  (6.14),  (6.15)  in  the  infinite  domain  by  the 
solution  of  equation  (6.14)  in  120  with  Robin  boundary  conditions  at  <9120.  We  have  established 
numerically  that  for  the  range  of  parameters  we  use  the  solution  of  the  problem  (6.14),  (6.15) 
is  close  in  12  to  the  solution  of  equation  (4.1)  in  the  bigger  rectangle  120  with  the  Robin 
boundary  conditions  at  its  sides. 

The  light  sources  are  located  in  several  positions  (xi,z)  =  (sj,10)  along  the  right-hand 
side  of  the  smaller  rectangle  12,  and  receivers  are  located  at  the  left-hand  side  of  12.  However, 
we  have  found  in  our  numerical  experiments  that  only  three  farthest  away  sources  are  typi¬ 
cally  useful  for  the  reconstruction.  We  use  three  (3)  sources  to  construct  an  approximation 
of  the  tail  functions  to  be  described  below.  Next,  we  use  two  (2)  sources  for  the  above  layer 
stripping  procedure  both  in  the  s-derivative  and  the  s-integral.  Although  it  is  possible  to 
use  more  data  or  light  source  positions,  but  our  numerical  experiments  showed  no  significant 
improvement.  We  have  also  introduced  the  multiplicative  random  noise  in  the  data. 

Following  the  convexihcation  method  described  above,  we  divide  the  domain  12  along  the  z 
axis  into  30  layers  within  the  interval  z  E  [^o,  zm  —  0\.  In  each  layer  we  then  approximate  the 
vector  function  p(z,  s )  along  the  z-axis  by  interpolations  using  quadratic  polynomials  (2.14). 
Coefficients  6j(s)  and  q(s)  are  determined  by  the  values  of  solutions  from  the  previous  layer. 
The  unknown  coefficient  aj(s)  is  determined  by  solving  equation  (2.33)  with  the  contractive 
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mapping  operator,  see  Theorem  2.1  in  section  2.  In  our  case  we  take  only  two  values  sq  and 
s2  =  s,  Si  <  S2  and  assume  that 


bi(s)  =  bi(s i),  Ci(s)  =  Ci(si),  a,i(s)  =  a*(s  1)  for  s  G  [si,  s2\  ■ 


In  x-direction  we  approximate  the  solution  using  Fourier  series.  Although  this  is  not 
completely  justifiable  (because  the  boundary  conditions  are  not  periodical),  it  does  provide 
a  reasonable  approximation  in  our  problem.  We  first  tried  the  number  65  of  terms  in  Fourier 
series  (2.21),  or  32  Fourier  modes  involving 


/iimx  7T  mx 

sm( — - — )  and  cos( — - — )  for  m  =  0,  ...,32. 


However,  this  case  was  unstable.  To  explain  the  latter,  we  note  that  the  problem  (2.22), 
(2.23)  inherited  the  instability  of  the  original  Cauchy  problem  for  the  nonlinear  integral 
differential  equation  (2.20a)  with  the  Cauchy  data  (2.20b).  It  is  well  known  that  the  Cauchy 
problem  for  an  elliptic  equation  is  unstable.  Hence,  in  our  case  the  numerical  error  will 
be  increasing  with  the  number  of  Fourier  modes.  In  order  to  reduce  the  numerical  noise 
generated  by  the  layer  stripping  reconstruction,  we  need  to  reduce  the  number  of  modes  in 
our  truncated  Fourier  series. 

Because  of  the  latter,  we  start  our  calculations  with  five  (5)  Fourier  modes  in  the  first 
layer  and  end  up  with  two  (2)  Fourier  modes  at  the  30st  layer  (the  right  edge).  To  do  this, 
we  change  initial  conditions  in  each  layer  as  follows 


b  intis') 


(i~  1)/A 

30  ) 


g  ~  l)h\ 

30  )  ’ 


where  bim{s )  is  the  component  of  the  vector  bi(s),  which  corresponds  to  either  of  functions 

.  7r  mi.  ,  nmx .  0 

sm( - )  or  cos( - ),m  =  3,4,5. 

5  5 


The  coefficient  Qm(s)  is  defined  similarly. 


6.4  Tails  v(x,  s)  in  (6.11) 

A  crucial  issue  in  our  problem  is  to  find  a  good  quality  approximation  of  the  function  v(x,  s) 
in  (6.11)  for  a  large  value  of  s.  In  the  case  of  inverse  problems  for  time  dependent  equations 
this  can  be  done  using  a  clear  asymptotic  behavior  of  the  Laplace  transform  of  the  solution 
of  the  forward  problem.  In  our  case,  however,  the  free  parameter  s  is  the  location  of  the  light 
source.  For  real  world  applications,  the  source  location  cannot  be  too  far  from  the  inclusion, 
both  due  to  the  restriction  in  size  and  the  limit  of  the  light  intensity.  We  have  undertaken 
a  substantial  effort  to  understand  the  behavior  of  solutions  when  locations  of  light  sources 
move  at  realistic  scales. 
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First,  we  consider  the  fundamental  solution  of  the  problem  (6.14),  (6.15)  for  the  case 
a(x,z)/D  =  k2.  This  solution  is 

u0  =  ^-K0(k\(x  s,  z  zm) |), 

ZTT 

where  Kq(z)  a  modified  Bessel  function.  It  is  well  known  that  the  asymptotic  behavior  of 
this  function  is 

„  ,  .  nr  ,, ,  ( _  „  /  1  \\  ,  ,  _ 


Kn(z)  = 


e~n  I  1  +  0  —  )),\z\  oo 


(6.17) 


Represent  now  solution  of  the  problem  (6.14),  (6.15)  with  the  condition  (6.16)  as  the  solution 
of  the  following  integral  equation 


u(x,z,s)  =  —K0(k\(x  -s,z-  zm) |) 

ZTT 

J  Ko  (kyj \x  -  02  +  (z  -  [a  (£,  v)  ~  k2]  u  (£,  V,  s)  d£dr>- 

n 

Let  S  ( x,z,s )  =  |  (x,z)  —  (s,  zm)|.  Introducing  the  function 

U (x,  z,  s )  =  2V2nSekS  ■  u(x,  z,  s ) 
and  taking  into  account  (6.17)  and  (6.18),  we  obtain  that 


e  kS  1  +  g{x,  z)  +  O  -  )  ,  S  ->  oo. 


(6.18) 


~  2^TkS  l  yK  ’  \S  J \  ’ 

The  function  g(x,z)  is  unknown  and  is  independent  of  S  as  S  — >  oo.  Hence,  we  obtain  for 
the  function  v  =  In  u 


(6.19) 


v(x,  z ,  s)  =  — kS  —  In  S'  +  g(x,  z)  +  O  I  —  J  ,  S'  — *  oo 


where  the  unknown  function  g(x,  z )  ^  g(x,  z ). 

We  approximate  the  function  g(x,  z)  by  two  different  methods  and  the  final  answer  is  the 
average  of  two.  We  start  at  z  =  zq  where  the  boundary  values  are  known.  We  decompose 
the  boundary  values  of  v  into 

5 

v(x,Zq,s)  =  Y  [&m(s)cos  +^m(s)sin  (~^~x)  >  (6.20a) 

m= 0 

Vz(x,  Zo,  s)  =  Y  [bizm(s)  cos  i~Yx)  +  bizl(s)  sin  {r?~x)  •  (6.20b) 

m= 0 
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Hence, 


bm(s)  =  ~^=  I  v(x,z0,s)cos  (-y-^)  dx, 


15 


(6.21a) 


bm\s)  —  y  j  v(x,  zq,  s)  sin  dx. 


15 


(6.21b) 


Therefore  it  follows  from  (6.19)  that  the  asymptotic  expansions  for  functions  bm\s)  are 


15 


b$(s)  =  —j=  j  -  [kS(x,z0,s)  +  \nS(x,  z0,s)]cos  dx  +  bi^ 


+ O 


S(x,Zq,s) 


The  function  g(x,  Zq)  can  be  approximated  by  a  truncated  Fourier  series  with  coefficient 
bm  (s)  and  bm\s)  and  similarly  for  gz(x,Zo).  To  approximate  numbers  bin  ,  we  take  three 
measurements  ranging  for  source  locations  Si,S2,Ss  that  are  far  enough  and  set 


£ 

k=\ 


15 


bm(sk)  +  y  y  (kS(x,z,sk)  +  lnS(x,z,sk))cos  (y-^)  dx 

5 


(6.22) 


where  numbers  bm\sk)  are  calculated  by  (6.21a,b).  We  do  similarly  for  other  coefficients  in 
(6.20a,b).  Note  that  in  (6.20a,b)  and  (6.22)  one  should  actually  put  sign  instead  of 
The  number  of  light  sources  N  —  3  is  taken  in  all  our  experiments  when  we  approximate  the 
function  g. 

However,  the  above  procedure  (6.19)-(6.22)  gives  us  the  value  of  the  tail  function  v  (x,  z,  s) 
and  its  derivative  vz  ( x,z,s )  in  (6.11)  only  at  z  :=  zq,  i.e. ,  v(x,zo,s) ,  vz  ( x,zq,s ) .  Equation 
(6.19)  provides  an  approximation  for  all  (x,  z)  £  O  if  we  simply  set  g(x,  z,  s )  =  g(x,  z0,  s ).  In 
our  numerical  experiments  we  found  that  this  is  insufficient.  Hence,  we  use  the  measurement 
data  from  a  different  angle,  which  enhances  our  numerical  results.  We  obtain  a  similar  tail 
function  using  the  measurement  data  at  the  lower  edge  of  H,  i.e.,  at  x  —  hem.  We  have 
decomposed  the  boundary  values  v  (5,  z,  s )  and  vz  (5,  z,  s )  into  a  Fourier  series  of  z  and  got 
a  second  tail  function  using  the  idea  similar  with  the  above.  Thus,  we  have  approximated 
v  (5,  z,  s) ,  vz  (5,  z,  s) .  Finally  we  set  for  the  tail 


v  (x,  z,  s )  :=  ^  [n  (x,  z0,  s)  +  v  (5,  z,  s)] , 

(6.23) 

(x,  Z,  s)  :=  ^  [vz  (x,  Zq,  s)  +  vz  (5,  z,  s)]  , 

(6.24) 
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In  our  numerical  experiments  we  have  used  the  value  of  A  =  50.  We  have  attempted  to  use 
the  tail  function  alone  from  (6.23),  (6.24)  for  the  reconstruction  of  a(x,  z).  The  tail  function 
has  provided  most  of  the  information  about  locations  of  the  inclusions.  These  locations 
were  reconstructed  precisely.  However  the  peak  value  of  the  reconstructed  coefficient  within 
inclusions  was  is  about  30%  off  the  original.  To  overcome  that,  we  follow  the  following  four 
steps  procedure: 

Step  1.  Use  the  tail  function  to  obtain  an  initial  guess  l-i(M(x,  z)  for  the  unknown  coefficient 

Step  2.  Re-scale  fia0(x,  z)  linearly  to  a  commonly  acceptable  biological  value.  In  our  case, 
we  set  / j,al(x,z )  =  /ia0(x,z)  if  fj,a0(x,z)  <  //,  below  the  threshold  and  /j,al(x,z)  =  a/j,a0(x,z) 
if  naQ(x,z)  >  where  ju  is  a  threshold  value.  The  scaling  constant  a  in  our  problem  is 
chosen  such  that  the  maximum  of  /xal(x,  z)  =  0.3  which  is  the  maximum  value  of  the  actual 
coefficient.  In  problems  where  the  actual  p,a{x,  z)  is  unknown,  we  use  biologically  well-known 
coefficient  for  that  specific  kind  of  medical  application.  The  value  of  Ji  ranges  from  0.14  to 
0.16  in  numerical  examples  we  calculated. 

Step  3.  Re-calculate  the  forward  problem  using  re-scaled  a±(x,z)  as  the  coefficient.  The 
function  v(x,z,s)  resulting  from  the  calculated  solution  is  used  then  as  the  tail  function 
instead  of  (6.23)-(6.24). 

Step  4-  After  we  have  obtained  the  function  u(x,z,s)  at  the  measurement  surface  at 
z  —  0  and  obtain  consequently  v  =  lnw  at  s  =  Si,  s2,  where  s2  >  Si,  we  set  s  :=  s2 ■  Then 
we  calculate  the  function  q(x,z,s i)  as 

v(x,z0,s2)-v(x,z0,s1) 
q{x,z  0,si)  = - . 

•‘>'2  —  -Si 

Similarly  for  the  derivative  qz(x,  z0,  Si).  This  way  we  obtain  initial  values  &o(si),co  (si)  •  In 
all  follow  up  calculations  on  all  layers  we  set  q(x,z,s )  :=  q(x,z,s i)  for  s  e  [si,s2]  •  I11  each 
layer  we  have  solved  equation  (2.33)  with  the  contractive  mapping  operator.  After  we  obtain 
q(x,  z,  si)  for  all  z,  we  update  v(x,  z,  si)  =  v(x,  z,  s2)  —  q(x,  z,  si)(s2  —  si) 

Step  5.  Reconstruct  the  coefficient  na(x,z)  from  the  function  v(x,z,s i)  using  a  finite 
element  code  of  [10]. 

Figures  15  and  16  display  some  results  of  reconstruction  from  the  transmitted  data  in 
the  2-D  case  with  x  =  (x,  z).  That  is,  the  source  was  running  along  the  line  {aq,  zm}  ,  zm  = 
const.  >  0,  x  G  [c,  d]  and  the  data  were  given  at  the  line  {x,  Zo}  ,  Zo  =  const.  <  zm,  x  G  [c,  d]  , 
as  well  as  at  the  line  x  =  {c,  z}  ,  z  G  [zq,  zm\  ,  where  the  data  at  the  second  line  were  used  to 
get  an  approximate  value  for  the  tail  function. 

6.5  Conclusions 

Although  the  case  of  the  running  source  is  not  covered  by  the  original  theory  of  the  convexi- 
hcation,  especially  the  issue  of  tails,  we  were  capable  to  adapt  the  convexikcation  scheme  to 
this  case.  Applications  include  imaging  of  land  mines  by  the  Electrical  Impedance  Tomog¬ 
raphy,  optical  imaging  of  targets  on  battlefields  covered  by  smogs  and  flames  and  medical 
optical  imaging. 
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Figure  15:  The  original  function  a(x,  z )  (top  left)  and  its  1-d  cross-section  (top  right)  via 
the  centerline  of  inclusions.  Other  1-d  cross-sections  of  Figs.  15  and  16  are  by  the  same  line. 
The  bottom  figures  display  the  reconstructiogiQfor  the  noiseless  data. 
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Figure  16:  Reconstruction  of  the  medium  of  Figure  15  (top)  with  2%  noise  in  the  data. 

7  Convexification  via  the  Cauchy- Riemann- like  Sys¬ 
tem  of  the  First  Order 


The  work  of  this  section  was  performed  by  Dr.  A.  Timonov,  who  was  working  during  2007 
under  a  subcontract  as  a  faculty  of  University  of  South  Carolina  Upstate.  He  has  proposed 
the  idea  of  replacing  the  second  order  PDE  (2.20a)  with  the  Cauchy-Riemann-like  system  of 
PDEs  of  the  first  order. 

Consider  again  the  Cauchy  problem  (2.20a,b).  To  simplify  the  presentation,  consider  the 
“tail  free”  case,  although  tails  were  incorporated  in  the  numerical  scheme.  We  rewrite  it  now 
in  the  2-Dimensional  case  as  follows 


S 

n 

S 

n 

qzz  +  Qxx  +  2  s2  qz 

-  qz  (x,  z,  r)  dr 

s 

+  2  s2qx 

-  /  qx  (x,  z,t)  dr 

s 

S 

p 

2 

S 

+2  s 

/  qx  (x,  z,t)  dr 

_  S 

+  2  s 

1 

h- 

In  addition,  the  following  Cauchy  data  are  given 


(7.1) 


q  (x,  0,  s)  =  ip 2  (x,  s ) ,  qz  (x,  0,  s)  =  ip2  (u  s) , 
for  (x,s)  £  (—A,  A)  x  (s0,s). 


(7.2) 
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Equation  (7.1)  contains  only  derivatives  of  the  function  q.  Hence,  introduce  two  new  func¬ 
tions  u  (. x ,  z,  s )  and  v  (x,  z,  s )  as 

u  (x,  z,  s )  :=  qz  (x,  z,  s ) ,  v  (. x ,  z,  s )  :=  (x,  z,  s ) .  (7.3) 


Then  (7.1)  leads  to 


uz  +  vx  +  2s2u 


s 

p 

S 

p 

—  /  u  (x,  z,  r)  dr 

J 

+  2  s2v 

—  /  v  (x,  z,  t )  dr 

J 

S 

S 

(7.4) 


S 

p 

2 

S 

p 

+2  s 

/  v  (x,  z,  t )  dr 

+  2s 

/  u  (x,  z1  t )  dr 

J 

J 

_  S 

_  S 

In  addition,  since  qxz  =  g23;,  then  (7.3)  implies  that 

vz-ux  =  0.  (7.5) 

And  by  (7.2) 

u(x,  0,  s)  =  -02  (x,  s) ,  v(x,  0,  s)  =  (p2x  (x,  s) .  (7.6) 

We  have  obtained  the  first  order  system  (7.4),  (7.5)  with  the  initial  conditions  (7.6).  If 
integrals  would  be  absent,  then  (7.4),  (7.5)  would  be  exactly  the  classic  Cauchy- Riemann 
equations  known  from  the  complex  analysis  for  real  and  imaginary  parts  of  an  analytic 
function. 

To  solve  the  problem  (7.4)-(7.6),  we  apply  an  analogue  of  the  convexihcation  method 
of  section  2.  However,  instead  of  using  quadratic  polynomials  on  each  layer,  we  use  linear 
functions  with  respect  to  z:  because  of  the  Erst  order  derivatives. 

Convergence  Theorem.  A  stability  and  convergence  theorem  similar  with  the  above 
Theorem  2.1  can  be  proven. 

It  is  natural  to  first  test  the  case  of  the  1st  order  equation  in  the  1-dimensional  case.  Three 
algorithms  were  tested  for  the  case  of  imaging  of  land  mines:  (a)  the  original  convexihcation 
algorithm  for  the  2nd  order  PDE,  (b)  the  above  new  version  of  the  convexihcation  for  the  1st 
order  PDE,  and  (c)  another  version  of  the  convexihcation,  which  is  based  on  the  introduction 
of  a  new  norm  using  the  Carleman  Weight  Function;  this  version  was  tested  for  the  1st  order 
PDE.  Figure  14  displays  some  of  obtained  numerical  results.  Since  the  performance  of 
versions  (a)  and  (b)  is  very  similar,  results  for  the  version  (a)  are  not  displayed  on  Figure 
17. 

Conclusion.  The  original  version  of  the  convexihcation  for  the  2nd  order  PDE  (section 
2)  and  two  new  versions  for  the  1st  order  PDE  provide  very  similar  numerical  results  in  the 
1-d  case.  Therefore,  the  activity  of  this  section  was  not  pursued  in  2008. 
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D  ry  sol 


Figure  17:  Numerical  comparison  of  the  performance  of  the  convexification  in  the  1-d  case 
for  the  first  order  equation.  Solid  line:  true  solution.  Bullets:  numerical  reconstruction  using 
the  first  order  PDEs  described  in  this  subsection.  Stars:  numerical  reconstruction  using  the 
re-normalization  and  the  first  order  PDEs.  Very  similar  results  were  obtained  by  the  original 
version  of  the  convexification  for  the  second  order  PDE. 
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8  Numerical  Solution  of  the  Inverse  Problem  of  Ther¬ 
moacoustic  Tomography 

Results  of  subsections  8.1  and  8.2  are  obtained  jointly  by  Dr.  C.  Clason  (Austria)  and  the 
PI  and  are  published  in  [10].  Results  of  subsection  8.3  are  obtained  jointly  by  the  graduate 
student  Mr.  A.  V.  Kuzhuget,  Drs.  S.  I.  Kabanikhin,  D.V.  Nechaev  (Russia)  and  the  PI, 
and  are  published  in  [11],  The  PI  is  the  thesis  advisor  of  Mr.  Kuzhuget  and  Kuzhuget  was 
partially  supported  by  this  grant. 

8.1  Theory 

Thermoacoustic  computed  tomography  is  a  new  imaging  method  that  uses  different  modal¬ 
ities  for  the  illumination  of  the  target  and  measurement  of  its  response.  The  target  is 
subjected  to  a  short  electromagnetic  impulse,  which  is  absorbed,  leading  to  a  temperature 
increase  and  hence  to  expansion.  This  induces  a  pressure  wave  in  the  target,  which  can 
be  measured  as  a  change  in  the  acoustic  held  outside  the  sample.  If  the  absorption  of  the 
electromagnetic  energy  is  spatially  varying,  the  resulting  wave  held  will  carry  the  signature 
of  that  inhomogeneity.  The  problem  is  hence  to  calculate  this  absorption  of  the  target  from 
time  dependent  acoustic  measurements  outside  it.  the  mathematical  statement  of  the  prob¬ 
lem  is  to  hnd  initial  condition  of  the  wave  equation  given  time  dependent  measurements  of 
both  its  solution  and  its  normal  derivative  at  a  surface. 

The  propagation  of  the  resulting  acoustic  pressure  held  u(x,t )  in  M3  is  governed  by  the 
equation 

-y-r Utt  =  Alt,  (x,  t)  e  M3  x  (0,  T)  (8.1) 

cz(x) 

with  the  initial  conditions 

u(x,  0)  =  fi(x),ut(x,  0)  =  f2  ( x ) ,  (8.2) 

where  the  function  f\  (x)  =  a  (. x )  k,  where  a(x)  is  the  spatially  varying  absorption  coefficient, 
which  is  unknown,  the  positive  constant  k  is  known,  so  as  the  spatially  varying  sound  speed 
c(x).  In  applications  to  thermoacustic  tomography  f2  (x)  =  0  and  only  the  function  fi(x)  is 
unknown.  However,  our  numerical  scheme  of  subsections  8.1  and  8.2  does  not  rely  on  any 
knowledge  of  functions  fi  (. x ) ,  f2  (x)  and  actually  determines  both  of  them  simultaneously. 
We  consider  the  following 

Inverse  Problem.  Let  SI  C  K3  be  a  bounded  domain  with  the  piecewise  smooth 
boundary  <90.  And  let  T  C  <90  be  a  part  of  the  boundary.  Denote  Tt  =  T  x  (0  ,T). 
Determine  functions  fi  (. x ) ,  f2  (x)  in  O  assuming  that  the  following  lateral  Cauchy  data 
(p0(x,t)  and  (p^Xjt)  are  known 

u  |rT=  (p0(x,t),dnu  |rT=  </h (x,t).  (8.3) 
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Denote  Qt  =  Ox(0  ,T)  ,St  =  <90  x  (0,  T ) .  Suppose  that  there  exists  a  function  g(x,  t )  G 
H 2  (Qt)  satisfying  boundary  conditions  (8.3),  i.e., 

9  |rT=^o 0M)A#  |rr=  ¥>i(aj, t).  (8.4) 

Denote  ^ 

w(t,  t)  =  u(x,  t)  -  #(t,  £),  L  =  -  A, 

ct(t) 

$0M)  =  ~L(g)(x,t). 

Then  (8.1)-(8.4)  imply  that 

Lw  =  <3>(x,  t)  in  Qy,  (8.5) 

iu  |rT=  dnw  |rT=  0.  (8.6) 

If  we  find  the  function  w  from  conditions  (8.5),  (8.6),  then  we  can  easily  find  unknown 
initial  conditions  fi(x)  =  w(x,  0)  +  g(x,  0),  /2(x)  =  w(x,  0)  +  g(x,  0).  Hence,  we  now  focus  on 
obtaining  a  good  approximation  of  the  function  w(x,t)  from  conditions  (8.5),  (8.6).  Clearly, 
the  problem  (8.5),  (8.6)  is  a  non-classical  one.  We  first  formulate  a  Lipschitz  stability 
estimate  for  this  problem,  which  is  Theorem  8.1.  This  theorem  is  proven  via  a  Carleman 
estimate  for  the  operator  L. 

Theorem  8.1  [10].  Assume  that 

Cmax  G"  C(x)  X'  Cmjn  X1  0,  Wx  G  D, 

there  exists  a  point  x0  G  R3  such  that 

^  +  (V  (<T2(t),t  >  0,  Vt  G  D. 


Let  r  =  maXj.go  \x  —  x0| .  Choose  a  number  /3  >  0  such  that 

\F&  <  /  _ a - - n/  -r  ,  Vt  G  D 

4(Cmin  +  r|Vc-2(T)|) 


and 

c2. 

/ 3  < 

3 

If  c=  const.  >  0,  then  these  conditions  reduce  to  c  >  f3.  Suppose  that  the  observation  time 
T  >  diameter  (fl)  /a fjd.  Then  there  exists  a  positive  constant  C  =  C  (Qt,  to,  cm;n,  cmax)  such 
that  the  following  Lipschitz  stability  estimate  holds 


HI hHQt)  -  C 


IMI 


L2  (Qt) 


UtII^1(5t)  +  1st  II L2(ST)  ’  G  R2  • 


This  theorem  guarantees  stability  of  the  numerical  method  proposed  below  for  the  case 
Tt  =  St-  To  solve  the  problem  (10.5),  (10.6)  numerically,  consider  the  Tikhonov  functional 
J£  (w)  , 

Je  H  =  | \\Lw  -  n2L2(QT)  +  |  \H\qr  ,  (8-7) 
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where  e  >  0  is  the  regularization  parameter  and 

IHIqr  :=  \\W\\l2{Qt)  +  Wwtt\\L2{QT)  +  IIAHIl2(qt)  • 

Denote  (,  )qR  the  scalar  product  corresponding  the  norm  ||.||qR  .  We  consider  this  functional 
over  the  space  H'l  (Qt)  ,  where 

Hq  ( Qt ,  hr)  =  {w  G  H2  ( Qt ,  IV)  :  w  |rT=  dnw  |rT=  0}  . 

By  the  variational  principle  the  minimizer  w£  of  the  functional  J£  ( w )  satisfies  the  following 
integral  identity 


LweLvdq  +  e  (w£,  v)=  QLvdq,  dq  :=  dxdt,  Wv  G  H q  ( QT ,  IV) 


(8.8) 


Qt 


Qt 


Hence,  we  need  to  solve  (8.8).  First,  we  formulate  the  existence,  stability  and  convergence 
theorems  for  (8.8).  Theorem  2  follows  from  the  classic  Ricsz  theorem. 

Theorem  8.2  [10].  For  any  function  G  L2  (Qt)  and  any  e  >  0  there  exists  unique 
solution  w£  G  Hq  (Qt,Tt)  of  the  problem  (10.8)  and 

„  C 

llw^ll^02(QT,rT)  <  II<|)IIl2(Qt)  ’ 


where  a  positive  constant  C  independent  on  £,  w£  and  <h. 

However,  Theorem  8.2  does  not  guarantee  existence  of  the  solution  of  the  original  problem 
(9.5),  (9.6).  To  formulate  convergence  and  stability  result,  we  assume  first,  that  T^  =  St- 
Also,  we  assume  by  the  Tikhonov  principle  [37]  that  there  exists  an  exact  solution  w*  G 
Hi  ( Qt,St )  of  the  problem  (8.5),  (8.6)  with  the  “ideal”  non-noisy  data  <h*.  In  addition,  we 
assume  that  the  actual  data  $  :=  <f>5  are  given  with  an  error,  i.e., 


'  L2(Qt) 


<*, 


where  <5  >  0  is  a  small  number  characterizing  the  level  of  the  error  in  the  data.  Let  wt  be 
the  solution  of  the  problem  (8.8)  with  <f>  :=  <f><5.  Theorem  8.3  guarantees  convergence  of  the 
regularized  solution  w&£  to  the  exact  solution  w*  as  £  :=  52  — >  0.  Theorem  8.3  follows  from 
Theorem  8.1. 

Theorem  8.3  [10].  Suppose  that  conditions  of  Theorem  8.1  hold  and  Tr  =  St-  Then 
there  exists  a  positive  constant  B  =  B  {Qt-  xn,  Cm\n,  cmaY)  such  that  the  following  estimate 
holds 


w£  -w 


I  Hi(QT) 


<  B 


52  +  £ 


1 1 2 

\\hhqt) 
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8.2  Numerical  Solution 


To  test  the  above  method  in  the  2-D  case,  we  took  =  (—3,3)  x  (—3,3)  and  T  =  7 
in  all  numerical  experiments  below.  We  have  used  the  Riesz-Galerkin  approximation  of  the 
function  w  (x,  y,  t )  via  cubic  5-splines,  which  provides  necessary  smoothness.  In  other  words, 
we  take  certain  partition  of  each  of  intervals 

rre  [-3,3 ],ye  [-3,3],te  [0,7] 
represent  the  function  w  (. x ,  y,  t )  as 


N  P 

w{x,y,t)  =  EE  Q'ijkBi  (x)Bj  (r/)5/j  (f ) 

i,j= 1  k= 1 

with  unknown  coefficients  a^.  As  test  functions  v  in  (8.8)  we  took  products  of  cubic  B- 
splines, 

buifci  (x,y,t)  =  Bix  (x)Bjl  (y)Bkl  (t). 

This  gave  us  a  linear  algebraic  system  Ma  —  F  with  respect  to  the  vector  a  =  {al3k)  ■  We 
have  solved  this  system  to  a  tolerance  of  10~6  by  a  stabilized  biconjugate  gradient  method 
(BICGSTAB,  provided  by  MATLAB).  Due  to  the  matrix  being  diagonally  dominant,  we 
employ  Jacobi  prescaling,  i.e.,  solving  DM  Da  =  DF  with  the  diagonal  matrix  Dtl  =  Mu  7  . 
This  has  proved  more  efficient  than  specialized  preconditioners,  since  the  matrix  is  diagonally 
dominant. 

To  test  the  stability  of  our  method  with  respect  to  the  random  noise  in  the  data,  we 
introduce  noise  in  calculated  spline  coefficients 

aijk  =  aijk  (l  +  aijki 

where  5  >  0  is  the  noise  level  and  £  is  the  random  variable  uniformly  distributed  over  [—1, 1] 
For  simulations  of  forward  problems  we  took  f2  (. x ,  y)  —  0  in  all  our  tests.  We  test  our 
algorithm  for  the  case  of  the  reconstruction  of  the  initial  condition  f\  (x,  y)  in  (8.2)  of  the 
form 

fi(x,  y)  =  sin(3x)  sin(3y)  exp  [-x2  -  y2]  , 
which  is  negligibly  small  outside  of  the  test  domain  D.  Denote 

fs(x,  y )  =  w£(x,  y,0)  +  g  (x,  y,  0) . 

Test  1.  Constant  coefficient  c.  We  take  Tt  =  St  and  c(x)  =  1.  Figure  18  shows  the 
slice  fe(x,  0)  for  £  =  10”3  for  various  noise  level  5  ranging  from  0  to  3.  One  can  see  that 
even  at  300%  noise  level  the  reconstruction  is  quite  good,  which  is  in  a  good  agreement  with 
the  previous  publication  [10].  The  Lipschitz  stability  estimate  of  Theorem  2  provides  an 
explanation  of  such  a  stability. 
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Figure  18:  Comparison  of  the  reconstructed  initial  condition  fie(x,y)  with  the  exact  one 
(solid  line)  in  a  homogeneous  medium  for  test  1  for  various  noise  levels  5  €  [0,3]  ,  i.e. ,  from 
0%  to  300%  noise.  The  regularization  parameter  5  =  1CT3.  Slices  of  f±£  (x,y)  are  shown. 


Test  2.  Non-constant  coefficient  c(x,y).  We  take 


1 

c2(x,y ) 


5 

2 


1 

12 


(*2  +  y2)  ■ 


(8.9) 


This  coefficient  satisfies  conditions  of  Theorem  8.1.  By  this  theorem  we  should  take  y/fi  < 
1/197,  which  leaves  us  with  T  =  1672.  However,  we  still  take  T  =  7,  since  we  believe  that 
the  estimate  of  Theorem  8.1  is  a  pessimistic  one.  Figure  19  displays  the  slice  f£(x,  0)  for 
e  =  10”3  for  various  noise  level  <5  ranging  from  0  to  3.  One  can  see  again  that  even  at 
300%  noise  level  the  reconstruction  is  quite  good.  Test  3.  Limited  boundary  data  with 
Ft  %  St-  Again  we  take  a  heterogeneous  medium  with  the  coefficient  c(x,y )  in  (8.9).  We 
take  T  =  14  and 

T  =  {(x,  y)  £  <90  :  {x  =  —3}  U  {y  =  —3}}  .  (8.10) 

Uniqueness  of  the  solution  can  still  be  proven.  Figure  20  shows  numerical  results  for  various 
noise  levels  5  £  [0, 1] . 


8.3  The  case  of  a  smaller  observation  time 

We  now  consider  the  same  Inverse  Problem  as  in  subsection  8.1  but  under  assumption  that 
at  least  one  initial  condition  at  t  =  0  is  known,  whereas  the  second  one  is  unknown.  We 
show  that  the  required  observation  time  T  can  be  lowered  by  the  factor  of  2  in  the  case 
when  in  (8.1)  c  =  1.  In  numerical  experiments  we  consider  the  case  of  unbounded  domain. 
Specifically,  quadrant.  It  was  shown  in  [11]  that  this  case  can  be  reduced,  under  certain 
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Figure  19:  Comparison  of  the  reconstructed  initial  condition  fie(x,y)  with  the  exact  one 
(solid  line)  in  a  homogeneous  medium  for  test  1  for  various  noise  levels  6  G  [0,1].  The 
regularization  parameter  6  =  10-3.  Slices  of  f\£  (x,  y) are  shown. 


Figure  20:  Results  for  the  case  of  incomplete  boundary  data  at  the  boundary  Y  in  (8.10). 
The  heterogeneous  medium  (8.9)  was  considered.  Slices  of  fi£  (x,y)  are  shown.  Due  to  the 
energy  loss  at  the  rest  of  the  boundary  <9D\r  the  resulting  solution  has  a  lower  amplitude. 
The  shape,  however,  is  preserved. 
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conditions,  to  the  situation  of  a  bounded  domain.  The  second  significant  difference  with 
the  above  is  that  instead  of  using  finite  elements  of  subsection  8.2,  we  use  in  our  numerical 
implementation  finite  differences.  This  enables  us  to  image  narrow  peaks,  which  model  5- 
function.  The  latter  was  impossible  when  using  finite  elements  of  the  previous  subsections, 
because  of  their  over-smoothness. 

Let  £2  C  Rn  be  a  convex  domain  with  a  piecewise  smooth  boundary  <9£2  and  2 R  be  the 
diameter  of  £2,2 R  =  rna xx,yeQ  \x  —  y\  ■  Let  T  =  const.  >  R.  Denote  Qt  =  £2  x  (0  ,T). 
Consider  the  elliptic  operator  L(x,t)  of  the  form 

n 

L(x,  t)u  =  A u  +  bj  (. x ,  t)  Uj  +  b0  (x,  t)ut  +  c  (x,  t)  u, 

3= i 

where  Uj  :=  dXju.  We  assume  that  all  coefficients  of  the  operator  L  belong  to  C  (Qt)  .  Let 
the  function  u  G  H 2  (Qt)  be  a  solution  of  the  hyperbolic  equation  in  the  cylinder  QT, 

utt  =  L(x,  t)u  +  F(x,t)  in  Qt,  (8.11) 

F  G  L2  (Qt)  with  initial  conditions 

u  (x,  0)  =  (p  (x) ,  ut  (x,  0)  =  ip  (x) ,  (p  G  H1  (£2) ,  t/j  G  L2  (£2) .  (8.12) 

In  this  subsection  we  consider 

Inverse  Problem  1.  Let  one  of  functions  p  or  vj  be  known  and  another  one  be  unknown. 
Determine  that  unknown  function  assuming  that  the  following  functions  /  and  g  are  given 

du 

U  |sT=  /  (x,  t),  —  | ST=  9  (x,  t ) ,  ST  =  <9D  x  (0,  T) ,  (8.13) 

where  v  is  the  unit  outward  normal  vector  at  <9£2.  We  call  the  problem  of  the  determination 
of  the  function  tp  the  “99— problem”  and  the  problem  of  the  determination  of  the  function  -0 
the  “^—problem”. 

It  follows  from  Theorem  8.1  that  in  the  case  T  >  2 R  one  should  not  assume  that  one 
of  functions  p  or  if}  is  known.  So,  we  are  interested  here  in  the  case  T  G  (R,  2 R) .  We  now 
specify  conditions  of  our  numerical  study.  Suppose  that  equation  (8.11)  is  homogeneous  with 
F  (x,  t)  =  0  and  it  is  satisfied  in  =  M2  x  (0,  T) .  Consider  the  quadrant  QU  =  {aq,  x2  >  0}  . 
And  also  consider  the  square  SQ  C  QU, 

SQ  (a)  =  {0  <  xi,x2  <  a}  . 

Suppose  that 

p(x)  —  ip  (x)  —  0  outside  of  SQ  (a)  .  (8-14) 

Then  the  energy  estimate  implies  that 

u  (x,  t)  =  0,  V  (x,  t)  G  {x  |  x  G  QU,  dist  (x,  SQ  (a))  >  T}  x  (0,  T) .  (8.15) 
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Figure  21:  Geometry  for  the  case  of  quadrant  (Inverse  Problem  2). 


Denote 

Tit  =  {x\  e  (0,  a  +  T) ,  x2  =  0}  x  (0,  T) , 

T2t  =  {%2  £  (0,  a  +  T) ,  X\  =  0}  x  (0,  T ) , 
r 3T  =  Id-'l  =  Cl  +  T,  £2  £  (0,  a  +  T)}  X  (0,  T) , 

r4T  =  {2^2  —  d  -\-  t,  x  1  g  (0,  a  +  t)}  x  (0,  r) , 
see  Figure  21.  Then  by  (8.15) 

Qu 

u  —  —  —  0  on  r3T  U  r4T.  (8.16) 

ov 

Hence,  we  focus  our  numerical  study  on 

Inverse  Problem  2.  Let  equation  (8.11)  be  satisfied  in  D/  with  initial  conditions  (8.12) 
satisfying  (8.14).  In  this  case  fl  :=  SQ  (a  +  T) .  Suppose  that  one  of  these  initial  conditions 
is  zero.  Determine  the  second  initial  condition,  assuming  that  functions  /  and  g  are  known, 
where 

du 

^  |riTur2T  f  (*^?  ^)  ?  ~Qy  lri^ur2T  9  (*^?  ^)  •  (8*17) 

It  was  proven  in  [11]  that  if 

T  >  (8.18) 

2  -  a/2 
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and  one  of  functions  p  or  ip  equals  zero,  then  the  following  Lipschitz  stability  estimate  is 
valid 

IMI h1(gt)  —  C  (ll/llifqrT)  IMIz^iy))  >  (8.19) 

where  rT  =  r1T  U  V2T  and  \\f\\m{rT)  =  \\f\\m(r1T)  +  ||/||Hi(r2T)  • 

We  consider  Inverse  Problem  1,  because  it  is  more  general  than  Inverse  Problem  2. 
Denote  Mu  =  utt  —  Lu.  To  solve  the  Inverse  Problem  1  numerically,  consider  the  Tikhonov 
regularizing  functional 


Je  (u)  —  II Mu  F\\L2[Qt)  +  £  ||w||i?2(QT) 


||  DPU  |  St  ~D,3f  ,  +  ||« 


2 

L2{St) 


v  \ St  dll 

2 


2 

L2  (St) 


(8.20) 


+X<p  \\MX,  °)  -  ^IIl2(Q)  +  Xy  IK®,  0)  -  v\\mm  ,Vu  e  H  (Qt) . 

Here  £  >  0  is  the  regularization  parameter, 

du 

Uu  I>St-—  I  St 

and  D 13 ,  \/3\  <  1  is  the  operator  of  (x,  t )  derivatives  with,  where  x- derivatives  are  those,  which 
are  taken  in  directions  orthogonal  to  the  normal  vector.  Also, 


Xtj) 


1  for  the  ip  —  problem 
0  for  the  p  —  problem 


i  Xm 


1  for  the  p  —  problem 
0  for  the  ip  —  problem 


Hence,  x^Xtp  —  0,  Xv  +  X*i>  —  1-  In  subsection  8.2  and  in  all  previous  works  on  the  QRM 
terms  in  the  second  line  of  (8.20)  were  absent  because  of  subtracting  off  boundary  conditions 
from  the  original  function  u.  Terms  in  the  third  line  of  (8.20)  were  absent  also,  and  they  are 
incorporated  now  to  emphasize  the  knowledge  of  one  of  initial  conditions. 

To  find  the  minimizer  of  Je  ( u ) ,  we  set  the  Frechet  derivative  of  this  functional  to  zero 
and  obtain  for  all  v  G  H 2  {Qt) 

J  MuMvdxdt  +  J  (p^vD^u  +  vuj  |,sT  dS  +  J  (vuuu)  \gT  dS 

Qt  St  St 


+Xt/j  /  [VuVu  +  uv }  (x,  0)  dx  +  x<p  /  ut{x,  0)vt{x,  0 )dx  +  £  [u,  v 


n 


FMvdxdt  +  £  (DPV  |sT)  D/3 fdS  +  +  {Vv  | St)  •  gdS 


Qt 


St 


<1 


St 


+X </,  J  [V^Vu(z,0)  +pv(x,0)\dx  +  xv  J  ipvt(x,0)dx. 


(8.21) 
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Riesz  theorem  and  (8.21)  imply 

Lemma  8.4.  For  any  vector  function  (F,  f,  g )  G  L2  (Qt)  x  Ft 1  (St)  x  L2  (St)  there  exists 
unique  solution  u£  G  Ft 2  (Qt)  of  the  problem  (3.2)  and 

\\uz\\h2{qt)  —  ~ffjf  (IIFIIl2(Qt)  ■*”  II/IIhrst)  +  \\9\\l2{st)  ”*■  Xy  IMIf/Rn)  +  W  ll^lli,2(n)J  • 

Setting  in  (8.21)  v  :=  u,  we  obtain  that  the  unique  minimizer  of  the  functional  J£  (u) 
satisfies  the  following  estimate 

1 1  ■ Mu  1 1 L  (Qt)  +  W  1 1  ■ u (x  >  °)  1 1  (n)  +  xv  1 1  (x ,  o)  1 1 12  (n) 

+  llM  IstIIhrst)  HWi/  IstIIl2(st)  (8.22) 

—  IIFIIl2(Qt)  +  II/II/TRSt)  +  IMI.L2(St)  W  IMIifqn)  +  Xt  llV’llijqn)  • 

To  prove  convergence  of  our  method,  we  need  to  derive  from  (8.22)  the  Lipschitz  stability 
estimate  for  the  function  u  in  the  Ft1  (Qt) -norm. 

Theorem  8.5  [11].  Let  C  Mn  be  a  convex  bounded  domain  with  the  piecewise  smooth 
boundary  and  let  T  >  R.  Suppose  that  the  function  u  G  Ft2  (Qt)  satisfies  the  inequality 

WMu\\ i/2 (Qt)  +  W  IKX >  °)  llifi(n)  +  II ut(x,  0)  || La(n) 

+  llM  I ST II J?1  (ST)  WUv  \stWl2(St)  —  ^r 
where  K  =  const.  >  0.  Then 

IMIhrQt)  +  Xv\W(x,  0)11^1(0)  +  Xylk  (®,0)||La(n)  <  CK. 

Theorem  8.5  enables  us  to  prove  convergence  of  our  method.  Following  the  Tikhonov 
concept  for  ill-posed  problems  [50],  we  first  introduce  an  “ideal”  exact  solution  of  either  p 
or  if  problem  without  an  error  in  the  data.  Next,  we  assume  the  existence  of  the  error 
in  the  boundary  data  /  and  g  and  prove  that  our  solution  tends  to  the  exact  one  as  the 
level  of  error  in  the  data  tends  to  zero.  We  consider  the  more  general  Inverse  Problem  1. 
Let  f*  G  Ft 1  (St)  and  g*  G  L2(St)  be  the  exact  boundary  data  (8.13),  F*  G  L2(Qt)  be 
the  exact  right  hand  side  of  equation  (8.11)  and  tp*  and  if*  be  exact  initial  conditions.  We 
assume  that  there  exists  an  exact  function  u*  G  Ft2  (Qt)  satisfying 

u*tt  =  L(x,t)u*  +  F*  (x,t)  in  Qt, 


with  initial  conditions 

u*  (x,  0)  =  ip*  (x)  ,  u*t  (x,  0)  =  ip*  (x) ,  ip*  G  Ft1  (12) ,  ip*  G  L2  (12) , 

du* 

u*  |sT=  f*  OM) ,  |sT=  9*  (x,t) , 
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where  p*  and  if*  are  exact  initial  conditions.  We  assume  that  the  real  boundary  data  in 
(8.13)  have  an  error,  so  as  the  given  initial  condition.  In  other  words,  we  assume  that 


11/  /  IIh^St)  +  Wd  9  II l2{st)  +  11-^  F  II L2(Qt) 

+W  \W  -  ^ll^i(n)  +  Xv  U  -  ^1l2(o)  < 
where  <5  >  0  is  a  small  number.  The  following  convergence  theorem  holds 

Theorem  8.6  [11].  Suppose  that  T  >  R.  Let  u£s  €  H 2  (Qt)  be  the  solution  of  the  QRM 
problem  (8.21),  which  is  guaranteed  by  Lemma  2.1.  Let  conditions  (5.1)- (5.4)  be  satisfied. 
Then  the  following  estimate  is  valid 


\\u  -  w*||hi(Qt)  +  W,  \\p  -  ^*||Hi(n)  +  W  U  -  V'l L2((i)  <  C  (5  +  y/e)  . 

In  our  numerical  study  we  have  considered  the  Inverse  Problem  2.  To  generate  the  data 
for  the  inverse  problem,  we  have  solved  the  Cauchy  problem 

ut.t  =  A u,  (x,  t)  G  M2  x  (0,  T) ,  (8.23) 

u(x,  0)  =  p  (x)  ,ut(x,  0)  =  if  (x) .  (8.24) 

In  our  numerical  experiments  if  (x)  =  0  for  the  p— problem,  and  p  (x)  =  0  for  the  ^—problem. 
Because  of  (8.15)  and  the  finite  speed  of  propagation,  we  use  in  our  solution  of  the  forward 
problem  zero  Dirichlet  boundary  condition  at  the  boundary  of  the  rectangle  (— T,  a  +  T)  x 
(— T,  a  +  T)  (Figure  18).  Hence,  we  solve  initial  boundary  value  problem  inside  of  this 
rectangle  for  equation  (8.23)  with  initial  conditions  (8.24)  and  zero  Dirichlet  boundary  con¬ 
dition.  In  all  our  calculations  we  took  a  =  1,  T  =  3.  The  square  SQ(a )  is  SQ(a )  =  SQ(1)  = 
(0, 1)  x  (0, 1) ,  the  domain  Q  is 

n  :=  (0,4)  x  (0,4) 

and  in  all  tests 

p(x)  —  if  (. x )  =  0  for  x  £  SQ(  1). 

We  have  solved  the  Cauchy  problem  (6.1),  (6.2)  via  finite  differences  using  the  uniform 
grid.  To  find  the  minimizer  of  the  functional  J£,  we  have  also  used  finite  differences  and  have 
minimized  the  resulting  functional  J£  with  respect  to  the  vector  {iq;mn}  ,  which  approximates 
values  of  the  function  u  at  grid  points.  Here  J£  means  the  functional  J£l  which  is  expressed 
via  the  finite  differences.  The  norms  \\uxi(x,  0)||L2(Q) ,  \\uX2(x,  0)||L2(^  in  \\u(x,  0)||^i^  in 

the  ^—problem  were  calculated  via  finite  differences.  As  to  the  term  \\D^u  |gT  s 

in  (8.20),  we  have  used  only  (3  —  0,  thus  ending  up  with  || u  |gT  —  /IIl2(5t)  c^scre^e 

sense) . 

To  minimize  the  functional  J£,  we  have  used  the  conjugate  gradient  method.  Derivatives 
with  respect  to  variables  Ukmn  where  calculated  in  closed  forms,  using  the  following  formula 


(dUkmn  _  t-  t-  t- 

~RfT  =  °kk°rnmVnni 
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where  5k k  is  the  Kronecker  symbol.  This  formula  can  be  conveniently  used  to  obtain  closed 
form  expressions  for  derivatives 

dJe  (u) 

dUkmn 


We  have  always  used  the  regularization  parameter  £  =  1CV6.  Larger  values  of  £  such 
as  1CT5  brought  lower  quality  results.  In  our  numerical  experiments  we  have  imaged  both 
smooth  slowly  varying  functions  and  the  finite  difference  analogue  of  the  6—  function.  Let 
(x\k,  X2r )  eObea  fixed  grid  point.  To  obtain  the  finite  difference  analogue  of  5  (aq  —  aq*,,  X2  —  X2r ), 
we  consider  the  following  grid  points  ( oqn ,  x2m)  and  model  the  function  5  ( oqn  —  aq*,,  a;2m  —  x2r) 
as 


&  (*^ln  •%! ki  J'2m  *^2 r) 


4hX2  hxi 


^nk^mrt 


where  the  multiplier  at  SnkSmr  is  chosen  such  that  the  volume  of  the  pyramid  based  on 
(xik-i,  X2r-i) ,  (xik-i,x2r+i) ,  (xik+i,x2r+i) ,  (xik+i,x2r-i)  equals  to  1.  Hence,  the  support 
of  the  function  S  ( aqn  —  Xik,X2m  ~  X2r )  is  limited  only  to  the  point  (aqn,aq m). 

Test.  The  (p— problem  with  two  6— functions.  We  now  consider  the  Inverse  Problem  2 
with  the  domain  as  with  'f(x)  =  0.  The  data  for  the  forward  problem  were  simulated  for 
the  case  of  two  S— functions 


p  (oq,  X2)  =  5  (x  1  —  0.4,  X2  —  0.4)  +  5  (x\  —  0.7,  X2  —  0.7)  (8.25) 

with  the  above  described  finite  difference  analogue  of  the  5—  function.  Figure  19  displays  the 
resulting  image  of  the  function  (8.25)  for  the  case  of  50%  of  the  noise  in  the  boundary  data, 
scatter  plot  mode  was  used,  squares  show  exact  height.  Very  similar  results  (not  shown) 
were  obtained  for  the  -0— problem  with  exactly  the  same  6—  functions  as  ones  in  (8.25). 

9  An  Inverse  Problem  for  a  Nonlinear  Parabolic  Equa¬ 
tion  with  Applications  in  Population  Dynamics  and 
Magnetics 

Results  of  this  section  are  obtained  jointly  with  Dr.  B.  Kaltenbacher,  University  of  Stuttgart, 
Germany  [14].  We  formulate  uniqueness  theorem,  outline  a  numerical  method  and  present 
some  numerical  results. 

Nonlinear  parabolic  PDEs  arise  in  a  large  variety  of  applications  ranging  from  combustion 
theory  via  environmental  pollution,  population  dynamics  and  nonlinear  magnetics  to  the 
theory  of  the  economic  growth.  Since  the  coefficient  depending  nonlincarly  on  the  solution 
of  the  PDE  is  often  not  accessible  to  direct  measurements,  its  determination  from  boundary 
measurements  is  an  important  task. 

We  consider  the  question  of  the  identifiability,  i.e.,  whether  this  coefficient  function  can 
be  uniquely  determined  from  the  given  data.  While  there  exist  results  on  the  situation  that 
the  coefficient  depends  on  values  of  the  solution  (as  it  is  relevant,  e.g.,  in  nonlinear  heat 
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Figure  22:  Test  of  subsection  8.3.  Exact  (red)  and  calculated  (black)  function  ip  with  50% 
noise  in  the  boundary  data.  The  function  %  =1  in  (8.20).  Scatter  plot  mode.  Squares  show 
heights.  Correct  heights  are  achieved. 
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conduction),  we  here  concentrate  on  the  problem  when  this  coefficient  is  a  function  of  the 
derivative  of  the  solution  (as  typical  for  nonlinear  magnetics,  (cf.,  e.g.,  [25])  or  where  the 
equation  is  in  non-divergence  form  (e.g.,  in  diffusion  models  for  population  dynamics,  cf. 

[24]). 

To  prove  our  uniqueness  theorem,  we  apply  the  above  mentioned  method  of  Carleman 
estimates.  The  majority  of  works  on  this  method  is  concerned  with  linear  equations  in  which 
unknown  coefficients  depend  on  spatial  variables.  Nonlinear  parabolic  equations  were  treated 
by  this  technique  in  [32],  [45],  [46]  and  Chapter  4  of  the  book  [41].  Our  inverse  problem  is 
a  problem  with  the  data  resulting  from  a  single  measurement  event.  In  the  case  of  multiple 
measurements  uniqueness  theorems  for  inverse  problems  for  nonlinear  parabolic  equations 
with  unknown  coefficients  depending  on  solutions  and  their  first  derivatives  were  proven  in 
[36].  The  case  of  single  measurement  data  for  such  equations  was  considered  in  e.g.,  [32], 
[41],  [45],  [46].  However,  in  these  works  the  unknown  coefficient  depends  on  the  solution  of 
the  original  parabolic  equation  (as  well  as  on  some  spatial  variables  in  the  multidimensional 
case).  The  case  of  the  dependence  on  derivatives  of  the  solution  was  not  considered.  Thus, 
it  is  an  important  new  element  here  that  the  unknown  coefficient  k  (■ ux )  is  involved  together 
with  its  first  derivative. 

9.1  Statement  of  the  Problem  and  Uniqueness  Theorem 

The  forward  problem  is 


ut  =  ( k  (ux)  ux)x  ,  (x,  t )  G  (0,  L )  x  (0,  T ) ,  (9.1) 

u  (x,  0)  =  r(x),  x  G  (0,  L) ,  (9.2) 

u  (0,  t)  =  f0  ( t ) ,  u(L,  t)  =  t  G  (0,  T) ,  (9.3) 

or  with  v(x,t)  =  ux(x,t),  i.e.,  u(x,t )  =  f*  v(£,t)  df  +  f0(t) 

vt  =  (k  (y)  v)xx  ,  ( x ,  t)  G  (0,  L)  x  (0,  T) ,  (9.4) 

v  (x,  0)  =  r'(x),  x  G  (0,  L) ,  (9.5) 

(k(v)v)x{0,t)  =  (k(v)v)x(L,t)  =  t  G  (0,  T) ,  (9.6) 

which  results  from  the  previous  setting  by  the  differentiation  with  respect  to  the  space 
variable.  Naturally,  it  is  assumed  that  the  function  k(z)  G  C 1  (M)  and 

k  (z)  >  k0  =  const.  >  0,  \/z  G  M.  (9.7) 


It  is  well  known  that  it  is  difficult  to  investigate  the  question  of  uniqueness  of  coefficient 
inverse  problems  for  parabolic  equations  unless  one  assumes  that  the  solution  of  the  forward 
problem  is  known  at  t  =  £  G  (0,  T).  In  the  latter  case,  however,  one  does  not  need  to  assume 
the  knowledge  of  the  initial  condition  at  (t  =  0} ,  see,  e.g.,  subsections  3.3.1  and  3.3.2  in 
[41].  In  other  words,  one  needs  to  assume  that  the  knowledge  of  the  initial  condition  (9.2) 
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is  replaced  with  the  knowledge  of  the  function  u(x,e)  for  an  arbitrary  £  G  (0,  T) ;  also  see 
below  for  a  physical  interpretation.  Because  of  this,  we  formulate  the  inverse  problem  as 
Inverse  Problem.  Assume  that  equation  (10.1)  is  satisfied  for  (x,t)  G  (0,  L)  x  (— c,  T) , 
where  c  G  (0,  T\  is  an  arbitrary  number,  the  function  u  (x,  0)  =  r{x)  in  (10.2)  is  known,  so 
as  boundary  conditions  /0  ( t )  and  f\  ( t )  for  t  G  (— c,  T ) .  Also,  assume  that  the  following  two 
functions  go(t )  and  g\{t)  are  given 

ux  (0,  t)  =  g0(t),ux  (. L,t )  =  gi(t),t  G  (-c,T) .  (9.8) 

Let  the  interval  (a,  b)  be  the  range  of  the  function  ux  (x,t)  for  (x,t)  G  (0,  L)  x  (— c,  T) ,  i.e., 
a  <  ux  (x,t)  <  b.  Determine  the  interval  ( a,b ),  as  well  as  the  unknown  coefficient  k  (z)  for 
z  G  (a,  b)  . 

The  following  uniqueness  theorem  holds. 

Theorem  9.1  [14].  Let  the  conditions  of  the  statement  of  the  Inverse  Problem  be  satisfied. 
Suppose  that  there  exist  two  functions  k\  (z)  and  (z)  satisfying  condition  (10.7)  and  such 
that 

zk\  (z)  +  hi  (z)  >  const.  >  0,  \/z  G  M,  i  —  1, 2. 

Let  Ui  (x,t)  and  u2  (x,t)  be  two  functions  satisfying  (10.1)  and  (10.3)  for  t  G  (— c,  T) ,  as  well 
as  the  same  Neumann  boundary  conditions  (10.8)  and  such  that  (see  (10.2))  iq  (x,  0 )=u2  (x,  0) 
r(x)  and 

DstDl.Ui  G  C  ([0,  L\  x  [— c,  T])  for  s  =  0,l;j  —  0, 1, 2,  3;  i  —  1, 2. 

Assume  also  that  in  (0 ,L)  x  (— c,  T) 

D2xUi  >  a  =  const.  >  0,  i  =  1, 2. 


In  addition,  let 


g'o(t)  >  0,t  G  (— c,  T) , 
gift)  =  const.  —  g,t  G  (— c,  T) 


and 

ki  (z)  =  k2(z)  ,z  G  (g0(-c),g0  (e)) , 

where  e  G  (0,  T)  is  an  arbitrary  number.  Denote  a  =  go(—c),b  =  g.  Then  the  range  of  both 
functions  u\x  and  u2x  coincides  with  the  interval  (a,  b)  and  k\  (z)  =  k2  (z)  for  z  G  (a,  b) . 


9.2  Application  Examples 

9.2.1  A  Diffusion  Model  in  Population  Dynamics 

In  this  section  we  discuss  the  assumptions  of  Theorem  9.1  for  an  example  from  population 
dynamics,  given  in  the  form  (10.4).  Aronson  in  [24]  has  proposed  the  diffusion  model 


vt  =  (4>(v)v)xx 


(9.9) 


in  the  absence  of  drift,  where  v(x,t)  is  the  population  density  at  location  x  and  time  t 
and  the  diffusion  coefficient  (j){y)  is  the  population  density  dependent  limit  (as  the  minimal 
individual  step  length  tends  to  zero)  of  the  second  moment  of  the  transition  probability 
between  different  locations.  The  function  0  plays  the  role  of  k  in  the  setting  of  Theorem  9.1. 
With  u(x,  t )  =  J*  v(£,  t)  d £  +  f0(t)  and  an  integration  with  respect  to  the  space  variable  this 
can  be  rewritten  as 

Ujt  (0(^x)^hr)x  C  (9.10) 

with  c  depending  on  time  only,  which  without  loss  of  generality  we  can  set  to  zero,  which 
yields 

x^)x  (9-H) 

so  that  with  k  =  0  we  get  the  PDE  considered  in  Theorem  9.1.  The  boundary  and  initial 
data  as  well  as  the  conditions  imposed  on  them  have  the  following  interpretation  in  the 
context  of  this  example: 

•  u(x,  0)  =  r(x),  i.e. ,  v(x,  0)  =  r'(x)  means  that  the  population  density  is  known  every¬ 
where  at  time  t  —  0,  e.g.,  from  a  census  at  that  instance  of  time. 

•  u(0,t)  =  f0(t )  is  satisfied  by  our  setting  of  the  integration  constants. 

•  u(L,t )  =  i.e.,  f()L  v(£,t)  =  fi(t)  —  fo(t)  means  that  the  total  population 

density  is  known  at  each  time  instance,  e.g.,  by  observation  of  the  in-  and  outflow  at 
the  boundary. 

•  ux(0,t)  =  u(0,  f)  =  g0(t)  with  g'  >  0:  The  population  density  at  location  x  =  0  is 
observed  for  all  time  instances  and  it  increases  with  time. 

•  ux(L,t )  =  v(L,t)  =  gi(t)  =  g  and  [go{e),g\  is  the  interval  on  which  k  can  be  identified 
according  to  Theorem  9.1:  The  population  density  at  location  x  =  L  is  observed  for 
all  time  instances  and  is  constant  in  time,  where  this  constant  should  be  as  large  as 
possible,  e.g.,  equal  to  some  saturation  value. 

•  uxx  =  vx  >  a  >  0:  The  population  density  increases  from  left  to  right. 

9.2.2  Nonlinear  Magnetics 

Quasistationary  magnetic  fields  can  be  described  by  a  subset  of  Maxwell’s  equations,  which 
leads  to  the  system  of  PDEs 


7u(  +  V  x  (iA7  x  u)  =  Jimp 

for  the  magnetic  vector  potential  u,  where  7  is  the  electric  conductivity,  /i  the  magnetic 
reluctivity,  Jimp  the  impressed  current  density  and  the  magnetic  flux  density  B  as  well  as 
the  magnetic  field  intensity  can  be  expressed  as 

B  =  V  x  u ,  H  =  uB  . 
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In  the  situation  of  high  magnetic  fields,  the  parameter  v  is  not  constant  but  depends  on  the 
magnetic  flux  density,  i.e., 

H  =  z/(B)B . 

Assuming  an  appropriate  geometry  (flat,  ring  shaped  probe  with  large  interior  radius)  we 
can  restrict  our  attention  to  the  spatially  one- dimensional  model  problem 

7 ut  -  (. v{ux)ux)x  =  0  , 

i.e.,  with  k  =  \v  we  arrive  at  the  PDE  considered  in  Theorem  9.1.  An  interpretation  of  the 
conditions  of  Theorem  9.1  in  the  context  of  this  example  can  be  made  as  follows: 

•  Neumann  boundary  data  ux(0,  t)  =  go(t),  ux(L,t )  =  gi(t)  here  have  the  meaning  of  a 
magnetic  flux  density,  which  can  be  extracted  from  measurements  of  the  magnetic  flux 
through  two  coils  positioned  on  both  sides  of  the  material  strip,  i.e.,  at  the  endpoints 
of  the  interval  [0,  L\.  The  magnetic  flux  can  be  controlled  via  the  impressed  current 
through  these  coils  such  that  it  is  monotonically  nondecreasing  (and  small)  at  the  left 
endpoint  as  well  as  constant  (and  large)  at  the  right  endpoint  to  achieve  g'0  >  0,  g[  =  0, 
and  a  wide  interval  (a,  b ). 

•  Via  the  relation  =  E  with  E  the  electric  field,  one  can  obtain  the  initial  and 
boundary  data  of  u  by  time  integration  of  electric  field  measurements. 

The  conventional  measurement  setup  for  determining  v  consists  of  just  an  excitation 
(and  measurement)  coil  wound  around  the  probe  and  does  not  enable  to  collect  all  the 
data  required  from  the  point  of  view  of  our  uniqueness  theorem  (especially  not  electric  field 
measurements  inside  the  probe  which  would  be  required  for  giving  initial  data  for  u) .  In 
this  sense,  there  is  a  gap  between  theory  and  practice.  Still,  the  model  problem  considered 
in  this  paper  gives  important  insight,  since  the  correct  form  of  the  PDE  with  a  coefficient 
depending  on  the  space  derivative  of  u  is  considered. 

9.3  Numerical  results 

In  this  subsection,  we  shortly  outline  a  numerical  method  for  determining  the  coefficient  k 
and  provide  some  computational  results. 

Here  we  make  use  of  a  special  time  discretization  based  on  trigonometric  polynomials, 
i.e.,  in  a  complex  valued  setting,  time  harmonic  functions  t  exp {int).  This  leads  to  a  so 
called  multiharmonic  formulation  of  the  forward  (cf.,  e.g.,  [25],  [35],  [37])  problem  as  follows: 
We  think  of  the  boundary  data  as  extended  periodically  and  sufficiently  smoothly  to  the 
larger  interval  [0,  T]  =  [0,  T  +  T\\  with  some  Tj  >  0,  and  make  an  ansatz 

N 

u(x,t )  ~uN(x,t)  :=  exp (mout)un(x) , 

n=—N 
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with  u  =  2tt/T.  (Note  that  continuation  of  the  boundary  values  beyond  T,  by  the  causality 
of  the  problem  does  not  influence  u  for  t  <  T .)  Inserting  into  the  PDE,  we  get 


N 


N 


Y  elJUJt  I  wjiij  -  /.'(  Y 


einuJtunx)ujx 


=  0 


j=—N 

Testing  this  with  the  functions 
and  using  the  orthonormality 


n=—N 


1 1— > ►  -=  e 

T 


—ihut 


_  /  g-tiart  vart  dt  =  S 

Tin 


where  Sij  is  the  Kronecker  Symbol,  yields 


~  ?  Jo  T.'L-N  e'(’  l|"‘  dt  =0 

V(e{-JV,:r.  ,,N}. 


(9.12) 


(9.13) 


To  eliminate  the  time  integral,  we  use  an  appropriate  approximation  for  k  that  enables  to 
take  advantage  of  the  orthogonality  (9.12).  For  this  purpose  we  make  use  of  a  polynomial 
ansatz 

p 

k(z )  ~  a pzp  , 

p= i 

which  is  justified  by  the  fact  that  in  the  applications  we  have  in  mind,  k  is  typically  a  smooth 
function.  Since  the  multinomial  theorem  yields 


with  the  multinomial  coefficients 


P 


pi 


iP—Nt  •  •  •  , Pn)  J  P-N]-  ■  ■  -Pn]-  ’ 
by  (9.12)  one  arrives  at  a  system  of  space  dependent  PDEs 


N  P 

uolui  —  (  EE  OLv(?i_ 

\n=—N  p= 0 


—r^nx 


=  0  l  e  {-N, 


(9.14) 


71 


where 


«=  E 


p  ex(p,s) 


T(p,s)  =  {P  =  (p~n,---,Pn)  C  INc 


Pn  =  P  A  UPn  =  S }  ’ 


i.e. ,  the  coefficients  cf_n  depend  nonlinearly  on  {u_Nx, . . .  ,unx )•  The  boundary  conditions 
and  measurements  become 


Mt*(0)  =  <?on  uix(L)  =  9n,  l  e  {-TV, . .  ■ ,  TV} 
«t( 0)  =  /oj,  ML)  =  fu,  le{-N,...,N} 


(9.15) 


fn  =  j,  J  exp (-iluT)fi(r)  dr, 

and  analogously  for  gt,  i  =  0,1.  Note  that  we  consider  the  Neumann  data  as  boundary 
conditions  (and  the  Dirichlet  data  as  measurement)  in  our  numerical  tests  below  in  order 
to  be  able  to  directly  prescribe  a  monotonically  increasing  g0  and  a  constant  g\ .  To  avoid 
singularity  of  the  system  due  to  the  elliptic  Neumann  problem  for  index  l  =  0,  we  make 
the  additional  normalization  assumption  u0  =  T  un(.,  t)  dt  =  0.  Physically,  this  corre¬ 
sponds  to  the  fact  that  in  experiments  only  higher  harmonics  (i.e.  multiples  nuj  with  |n|  >  1 
of  the  basic  frequency  u)  appear.  From  a  mathematical  point  of  view  this  normalization 
corresponds  to  a  T  periodicity  assumption  on  the  antiderivative  of  uN ,  while  the  multihar¬ 
monic  ansatz  implies  periodicity  of  uN  itself.  To  keep  compatibility,  we  add  this  assumption 
of  vanishing  time  integral  also  to  the  conditions  imposed  on  the  periodic  extension  of  the 
boundary  values  ./o,  /i,.9o,.9i- 

Since  k  is  real-valued,  we  expect  to  get  a  real-valued  solution  and  therefore  additionally 
stipulate 

u-n  =  un.  (9.16) 

The  inverse  problem  of  reconstructing  k  from  boundary  measurements  of  u  in  this  kind  of 
discretization  amounts  to  determining  the  polynomial  coefficients  ao, . . . ,  ap  from  boundary 
data  of  ?h,  l  —  0, . . . ,  N.  Writing  this  as  a  system  of  equations 

K(a)  =  V 

where  :  a  ^  (unJ0),unJL))n^)N  with  (u„)„=0  solving  (9.13),  (9.15),  (9.16),  we  can 
apply  Newton’s  method 

am+1  =  am  +  9(3  with  F£\am)P  =  V~  F^{am) , 

to  iteratively  recover  a.  In  here  6  is  an  appropriately  chosen  stepsize  that  guarantees  strictly 
monotone  decrease  of  the  residual.  Keeping  in  mind  the  fact  that  the  original  infinite  di¬ 
mensional  inverse  problem  of  recovering  k  is  ill-posed  in  the  sense  of  instability,  we  here  rely 
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Figure  23:  A  typical  B-H  curve  for  iron. 


on  the  regularizing  effect  of  discretization.  More  precisely,  we  stabilize  the  problem  by  pro¬ 
jection  in  data  space,  since  our  special  time  discretization  corresponds  to  an  L2  projection 
in  data  space  to  the  subspace  spanned  by  the  first  2N  +  1  trigonometric  polynomials.  These 
L2  projections  converge  pointwise  to  the  identity  due  to  the  fact  that  {t  i— >•  -U \  j  6  7L} 
forms  a  complete  orthonormal  system  in  L|.(0,  T).  To  fully  discretize  this  method,  it  remains 
to  define  a  space  discretization  of  the  space  dependent  functions  ui,  which  we  do  by  using 
finite  elements  on  a  uniform  grid. 

The  application  we  have  in  mind  here  is  nonlinear  magnetics,  where  the  nonlinearity  is 
usually  considered  in  terms  of  the  B-H  curves  a  where  a~l(z)  =  k(z)z.  B-H  curve  means 
a  plot  of  the  magnetic  flux  density  B  over  the  magnetic  field  intensity  H.  A  typical  B-H 
curve  for  iron  is  displayed  in  Figure  23,  where  the  ratio  between  the  slope  at  B  =  0  and  the 
slope  at  B  —  Bsat  amounts  to  the  relative  reluctivity  /irel  ~  ^A_..  por  testing  the  proposed 
method,  we  constructed  two  test  problems  with  similar  behavior.  Upon  renormalization 
to  =  1,  the  values  H™[m  =  3.6  and  H™°fm  =  2.4  as  appearing  in  the 

synthetic  test  example  in  Figure  24  indeed  turn  out  to  be  in  a  realistic  range. 

For  the  test  example  of  Figure  20,  Figure  21  shows  the  iteration  history  of  the  proposed 
method  for  exact  and  for  randomly  perturbed  data  with  a  noise  level  of  one  per  cent,  as 
compared  to  the  true  curve  aex  in  black.  Note  that  in  our  numerical  tests,  the  starting  value 
for  k  is  just  a  constant  function,  corresponding  to  the  fact  that  in  practice  the  constant 
reluctivity  coefficient  for  small  magnetic  fields  is  typically  known. 
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Figure  24:  Starting  value  a0  and  Newton  iterations  a\ ,  a2,  a3,  a4  for  a,  where  a  1  (^)  =  ^/c(^). 
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10  The  First  Result  on  the  Second  Generation  of  Glob¬ 
ally  Convergent  Numerical  Methods 

This  activity  is  an  ongoing  joint  effort  of  Dr.  L.  Beilina  (Norway)  and  the  PI  [12].  Prior  to  this 
activity  Dr.  Beilina  has  developed  (with  co. -authors)  the  so-called  adaptivity  technique  for 
inverse  problems  [26-28] .  The  idea  is  to  develop  the  second  generation  of  globally  convergent 
numerical  methods.  The  first  step  of  such  methods  is  the  same  as  in  the  convexihcation:  the 
transformation  procedure  which  leads  to  the  nonlinear  integral  differential  equation  (2.20a). 
However,  the  question  of  numerical  solution  of  this  equation  is  the  most  challenging  one 
in  the  entire  approach.  Thus, the  main  difference  with  the  convexihcation  is  in  the  solution 
of  this  equation.  Instead  of  using  the  layer  stripping  procedure  with  respect  to  a  spatial 
variable,  we  now  use  the  layer  stripping  with  respect  to  the  pseudo  frequency  s  via  going 
from  the  truncation  high  pseudo  frequency  in  the  downward  direction.  On  each  small  s-layer 
we  solve  the  Dirichlct  boundary  value  problem  for  an  elliptic  equation  via  the  FEM.  Since 
the  derivative  with  respect  to  s  is  not  a  part  of  the  differential  operator  (unlike  derivatives 
with  respect  to  spatial  variables  of  the  convexihcation),  we  expect  a  better  stability  of  this 
second  generation  methods.  Further,  the  Dirichlet  boundary  value  problem  presumes  that 
the  “measurements”  are  taken  over  the  entire  boundary.  At  the  same  time,  in  problems 
of  the  interest  to  the  Army  only  back  rehected  data  are  given.  Nevertheless,  it  seems  that 
this  new  approach  can  be  extended  to  this  case.  To  do  this,  one  should  impose  radiation¬ 
like  boundary  conditions  (i.e.,  Robin  conditions  actually)  on  the  rest  of  the  boundary.  In 
this  case  one  should  replace  the  Dirichlet  boundary  problems  with  “mixed”  boundary  value 
problems.  In  other  words,  Dirichlet  conditions  will  be  assigned  on  the  “accessible”  part  of 
the  boundary  and  Robin  conditions  on  the  rest.  However,  the  latter  development  is  a  matter 
of  future  effort,  since  only  the  first  result  in  this  direction  is  obtained  so  far  [12]. 

Another  essentially  new  and  important  element  here  is  the  presence  of  s-dependent  Car- 
leman  Weight  Functions  (CWFs)  in  the  numerical  scheme,  s-dependent  CWFs  are  present  in 
our  numerical  scheme.  Until  now  CWFs  have  been  widely  used  for  proofs  of  unique  continu¬ 
ation  and  conditional  stability  results  for  ill-posed  Cauchy  problems  for  PDEs,  as  well  as  for 
multidimensional  CIPs  with  the  single  measurement  data,  see,  e.g.  [29,30,32,38,39,41,45].  In 
this  capacity  CWFs  were  dependent  on  spatial  variables,  since  they  have  provided  weighted 
estimates  for  differential  operators.  In  the  convexihcation  CWFs  also  depend  on  the  spatial 
variable  z.  However,  CWFs  of  this  new  method  are  used  for  integral  Volterra-like  operators, 
they  are  involved  in  the  numerical  scheme  and  depend  on  the  pseudo  frequency  s,  rather 
than  on  a  spatial  variable. 

We  now  state  the  inverse  problem  a  little  bit  differently  than  in  section  2.  Namely,  we 
use  only  one  boundary  condition  at  the  entire  boundary.  And  also  we  do  not  require  that 
the  domain  O  be  a  rectangular  prism. 

Inverse  Problem.  Let  D  C  R3  be  convex  bounded  domain  with  the  boundary  G  C2. 
Suppose  that  one  of  coefficients  of  the  equation  (2.9a)  is  unknown  in  f2  and  known  in  M3\f! 
and  all  other  coefficients  are  known  everywhere.  Determine  that  unknown  coefficient  for 
i6fl,  assuming  that  the  following  two  function  tp  (x,  s )  is  known  for  a  single  source  position 
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xq  n 


(10.1) 


w  (x, ,  s)  =  (p  (x,  s ) ,  for  (x,  s )  G  — <912  x  [so,  s]  , 

where  so  and  s  are  certain  positive  numbers. 

Denote 

d  fhup  (x,  s) 

10.1  Layer  Stripping  With  Respect  to  the  Pseudo  Frequency 

We  approximate  the  function  q  (x,  s )  in  (2.20a)  as  a  piecewise  constant  function  with  respect 
to  the  pseudo  frequency  s.  That  is,  we  assume  that  there  exists  a  partition  s  =  S]y<  sjv_i  < 
...  <  Si  <  s0  =  s  of  the  interval  [s,  s]  with  the  sufficiently  small  grid  step  size  h  =  Sj_i  —  st 
such  that  q  (x,  s )  =  qn  (x)  for  s  G  (sn,  sn-i]  •  Hence, 

^  re— 1 

/  X7q  (x,  t )  dr  =  (sn_i  -  s)  Vqn  (x)  +  h  ^  Vqj  (x)  ,  s  G  (sn,  s„_i]  . 

,  i=1 


We  approximate  the  boundary  condition  (10.1)  as  a  piecewise  constant  function, 

qn  (x)  =  i\)n  (x)  ,  x  G  <912, 

where 


(10.2) 


Sn—  1 


=  ^  J  (x,  s)  ds. 

Sn 

Hence,  for  .s  G  (sn,  sn_i]  the  equation  (2.20a)  can  be  rewritten  as 

(re-1 

h^2^qj  (x)  ) 

j=i 

+2  (s2  -  2s  (sn_i  -  s))  V<?„ W  (x) 

=  2  (sn_!  -  s)  [s2  -  s  (sn_i  -  s)j  iy  qn)2  ~  2sh2  V<^ 


'  n—  1 


[X 


(10.3) 


o=i 


n— 1 


+4sW  (x)  Vgj  (x)j  -  2s  |W  (x)\ 2 ,5  G  (5n_i,5n] 

The  equation  (10.3)  is  nonlinear,  and  it  depends  on  the  parameter  s,  whereas  the  function 
qn  (x)  is  independent  on  s.  This  discrepancy  is  due  to  the  approximation  of  the  function 
q  (x,  s )  by  a  piecewise  constant  function.  Although  it  seems  that  the  equation  (10.3)  is  over- 
determined  because  the  function  qn  (x)  is  not  changing  with  the  change  of  s,  variations  of 
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s-dependent  coefficients  of  (10.3)  are  small  over  s  G  [sn,  sn_i) ,  because  this  interval  is  small. 
This  discrepancy  is  actually  helpful  for  our  method,  because  it  enables  us  to  “mitigate”  the 
influence  of  the  nonlinear  term  (Vgn)2  in  (10.3)  via  introducing  the  s- dependent  CWF. 

In  addition,  we  add  the  term  —  eqn  to  the  left  hand  side  of  equation  (10.3),  where  £  >  0  is 
a  small  parameter.  We  are  doing  so  because,  by  the  maximum  principle,  if  a  function  p(x,  s) 
is  the  classical  solution  of  the  Dirichlet  boundary  value  problem 

Ln(p)  -ep  =  f(x,s )  in  Q,p  \aa=Pb(x,s), 


then 


max  \p\  <  max  max  \pb\  ,  £  1  max  |/| 
o  9n  ci 


,  Vs  G  (sn—i,sn  . 


(10.4) 


On  the  other  hand,  if  e  =  0,  then  an  analogue  of  the  estimate  (10.4)  would  be  worse  because 
of  the  involvement  of  some  other  constants.  Therefore,  it  is  anticipated  that  the  introduction 
of  the  term  —eqn  should  provide  a  better  stability  of  our  process,  and  we  indeed  observe  this 
in  our  computations. 

Introduce  the  s-dependent  CWF  Cn, \  (s)  by 


f^n,A  (®)  exp  A  |s  Sn__i|  ,s  G  (sn, Sn_i  , 


(10.5) 


where  A  >>  1  is  a  parameter.  Multiply  both  sides  of  (10.3)  by  this  CWF  and  integrate  with 
respect  to  s  over  the  interval  [sn,  sn_i] .  We  obtain 

+  AinVqnW  —  sqn 

V qt  (x)^  (10.6) 


+2A2jUVV 

where 


n—  1 


]TVg,:  —  A2,n  (VC)2  ,  n  =  1, ...,  N, 


i=  1 


n—  1 


Ln  (qn)  ■=  A qn  -  Ahn  h  ^  Vg,  Vg,, 


i= 1 


2-yA  (Vg„)2  -  A2;nh2  (  ^ 


'  n—1 


„  i=  1 


I0:=I0{\,h) 


Sn  —  1 

J  Cn,x  (s)  ds 

Sn 


1  -  e~Xh 
A 


d  1  ,n  ■  d\  n  (A,  ^l) 


Sn—  1 

J  (5)1—1  5)  [s  s  (sn_i  s)J  Cn  \  (s)  ds , 

Sn 


^4l,n  •  Ai,n 


Sn—  1 


Sn 


2s  (sn—  i  <§))  Cn^\  (^)  ds^ 
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^4-2 ,n  •  -'4.2, n  (A, 


2 

A 


•Sn  —  1 

J  sCn.x  (s)  C?S. 

Sn 


Thus,  we  have  obtained  a  boundary  value  problem  (10.2),  (10.6)  for  a  coupled  system  of 
nonlinear  elliptic  PDEs  with  respect  to  the  vector  function  (g1; ...,  q^).  In  this  system  the 
tail  function  V  is  also  unknown  and  should  be  treated  differently.  An  important  observation 


for  A h,s  >  1.  Therefore  by  taking  A  >>  1,  we  mitigate  the  influence  of  the  nonlinear  term 
with  (Vg„)2  in  (10.6). 

In  computations  the  above  integrals  with  the  CWF  should  be  calculated  in  closed  forms. 
This  is  because  the  function  Cn_\  (s)  changes  rapidly  for  large  A,  which  means  that  the 
integration  step  size  should  be  taken  too  small.  In  principle,  one  can  decrease  the  step 
size  h  in  the  s-direction  instead  of  using  the  CWF.  However,  the  introduction  of  the  CWF 
provides  more  flexibility  for  the  choice  of  parameters  for  computations,  since  parameters  h 
and  A  are  independent,  as  long  as  Xh  >  1.  In  addition,  taking  h  too  small  would  increase 
the  computational  time,  because  one  would  need  to  compute  too  many  functions  qn  then. 


10.2  The  Algorithm 

The  above  considerations  lead  to  the  algorithm  described  in  this  section.  Below  Ck+a  (O) 
are  standard  Holder  spaces,  where  k  >  0  is  an  integer  and  a  G  (0,1).  Denote  \f\k+a  = 
||/||<7k+a(n)  ,  V  /  G  Ck+a  (O)  .  Our  algorithm  reconstructs  iterative  approximations  cn^  (. x )  G 

C°  (O)  of  the  function  c  (x)  only  inside  the  domain  O.  On  the  other  hand,  to  iterate  with 
respect  to  tails,  we  need  to  solve  the  forward  problem  (2.9a,b)  in  the  entire  space  M3.  To  do 
this,  we  need  to  extend  each  function  cn^  (x)  outside  of  the  domain  O  in  such  a  way  that 
the  resulting  function  cnjfc  G  Ca  (M3) ,  cn^  >  d\  in  0  and  cn^  =  2di  outside  of  O.  So,  we 
first  describe  the  procedure  of  such  an  extension  and  this  procedure  is  a  rather  standard  one. 
Choose  a  smaller  convex  subdomain  IT  C  O.  Choose  a  function  y  (x)  G  C 1  (M3)  such  that 

(  1  in  O', 

X  (x)  —  s  between  0  and  1  in  0\0', 

[  0  outside  of  O 

The  existence  of  such  functions  y  (x)  is  well  known  from  the  Real  Analysis  course.  Define 
the  target  extension  of  the  function  cnj.  as  (x)  :=  2d\  (1  —  y  (x))+y  (x)  cn ^  (x) ,  \/x  G  M3. 
Hence,  and  cn^  (x)  =  const.  =  2 d\  outside  of  the  domain  O  and  'cn^  G  Ca  (M3).  Furthermore, 
if  cn,k  (x)  >  di  in  O,  then  also  cn,k  (x)  >  d\  in  M3.  Indeed, 

Cn.fc  ( x )  -  di  =  (1  -  y  (x))  d i  +  y  (x)  ( c  (x)  -  d i)  >  0. 
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Step  l1.  Choose  an  initial  tail  function  V\A  (. x )  G  C2+a  (hi).  Choose  a  large  parameter 
A  >>  1  and  a  small  parameter  £  G  (0, 1) .  To  compute  the  first  approximation  qi ;i  for  the 
function  q\  with  this  tail,  solve  the  following  Dirichlet  boundary  value  problem 

Aqi, i  +  AltlVq1AVV1A  —  eqiA  =  —A2A  (Wqi)2  , 
qi,i  =  ip i  (x) ,  x  G  <9fL 

By  the  classic  Schauder  theorem  this  problem  has  unique  solution  q\,\  G  C2+a  (fl)  .  Re¬ 
construct  an  approximation  C\A{x)  G  Ca  (Q)  for  the  unknown  coefficient  c(x)  using  the 
function  g1;1  (x)  and  formulas  (2.13),  (2.18),  where  x  G  fi,  V  (x,s)  :=  :=  Si. 

Next,  construct  the  function  c.\A  G  Ca  (M3) . 

Step  lfc,  k  >  2.  Suppose  that  the  function  ci^-i  G  Ca  (O)  is  reconstructed,  Ci^-i  (x)  > 
d\  in  and  the  function  G  Ca  (M3)  is  also  constructed.  Solve  the  forward  problem 

(2.6),  (2.7),  where  c(x)  :=  Ci,k-i  (x),  s  —  s.  Let  W\ ,k(x,s)  be  the  solution  of  this  forward. 
Update  the  tail  function  as 

Ui.fe  (x)  =  ~2  hr  wltk(x,s)  G  C2+a  (Q)  . 

Next,  solve  the  boundary  value  problem  for  the  equation 

Aqitk  +  ^4i,iVq,i,fcVUi,fc  —  £qi,k  =  2  — —  (V^i^-i)  —  A2A  (VVi^) 

to 

with  the  boundary  condition  (10.2)  (at  n  —  1).  We  obtain  the  function  gljfc  G  C2+a  (O)  . 
Reconstruct  a  new  approximation  C\A  G  Ca  (fl)  for  the  unknown  coefficient  similarly  with 
the  Step  l1  in  which  U1;1  (x,s)  is  replaced  with  V\A  (x,s).  In  addition,  construct  the  func¬ 
tion  ci,fc_i  G  Ua  (M3) .  Make  several  steps  l1,  l2, lmi.  As  a  result,  we  obtain  functions 
9i,Ci,Ci)mi,Vi)mi  where 

<h  '■=  9i,mi  G  C2+q  (U)  ,  Ci  :=  Cijmi  G  CQ  (hi)  ,  Ci)TO1  G  Ca  (M3)  ,  Ui,mi  G  C"+°  (fi)  . 

Step  n1.  Having  functions  qi,..,,qn-i  G  C2+a  (fl)  and  the  tail  function  G 

C2+a  (fi)  ,  set  gn,0  (x)  :=  gn_i  ( x )  in  H, 

hn,l  (*^)  •  Ki-l,mn_i  (*^)  ■>  x  G  12. 


Solve  the  following  Dirichlet  boundary  value  problem  for  the  function  qnA  (in  (10.8),  (10.9) 
the  vector  function  (qn_k,  qn,o,  V^*,)  should  be  replaced  with  (qn,i,  qn,k-i,  Ki,i)) 


A 


£qn,i  T  Ai)nVg„,fc  •  W„)fc 


2 -A-  (Vgn,fe_i)2  —  A2^nh2 

to 


(10.8) 
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A2 ,n  (V^fc)2  , 


(re- 1 

h  (x) 

3= 1 

9re,l  (®)  =  V’n  (X)  ,  ®  £  <9ft.  (10.9) 

Hence,  we  obtain  the  function  gn  l  G  (72+a  (ft)  .  Reconstruct  a  new  approximation  cnji  G 
C'Q  (ft)  for  the  unknown  coefficient  c(a:)  using  the  function  qn>i(x),  functions  qi,...,qn-i 
and  formulas  (2.13),  (2.18),  where  x  €  ft,  V  (. x )  :=  Vn>i  (x),s  :=  sn.  Next,  construct  the 
function  cn,i  E  Ca  (M3) . 

Step  nk,  k  >  2.  Suppose  that  the  function  cn)k- 1  G  Ca  (ft)  is  reconstructed,  cn^-i  (x)  > 
di  in  ft  and  the  function  'cn^-i  G  Ca  (M3)  is  also  constructed.  Solve  the  forward  prob¬ 
lem  (2.9a,b),  in  which  c(x)  :=  'cn,k- l  (x)i  s  —  s.  Let  wn^{x,s)  be  the  solution  of  this  for¬ 
ward  problem.  Llpdate  the  tail  function  as  Vn^(x)  =  s-2  In  wUjk{x,  s)  E  C2+a  (ft)  .  Next, 
solve  the  boundary  value  problem  (10.8),  (10.9).  Reconstruct  a  new  approximation  cn^  E 
C°  (ft)  for  the  unknown  coefficient  similarly  with  the  Step  n 1  in  which  14, i  ( x ,  s)  is  re¬ 
placed  with  14, fc  (x,s).  Make  several  steps  n1,^2,  ..,nmn.  As  a  result,  we  obtain  the  functions 
<?n,cn,cn,K,mn,  where 


qn  ■=  qn,mn  G  <72+Q  (ft)  ,  cn  :=  cn,mn  G  <7C 


G  C0 


n,mn 


X 


E  C2+c 


If  functions  cn(x)  did  not  yet  converge,  then  proceed  with  Step  (n  +  1)1,  provided  that 
n  <  N .  However,  if  either  functions  cn(x)  converged,  or  n  =  N ,  then  stop.  The  convergence 
criterion  for  functions  cn(x),  which  we  have  established  in  our  computational  experiments,  is 
described  in  sub-subsection  10.4.2.  In  principle,  however,  there  might  be  several  convergence 
criteria,  which  indicates  that  a  better  one  might  be  found.  Thus,  we  do  not  specify  such  a 
criterion  here. 


10.3  Global  convergence  theorem 


First,  we  reformulate  the  Schauder  theorem  in  a  way,  which  is  convenient  for  our  case. 
Assuming  that 

s>l,Xh>l  (10.10) 

and  (10.7)  and  formulas  for  numbers  Aitn,A2>n,  we  obtain  that  |Aijn|  <  6s2  and  \A2)U\  <  2s 
for  all  n  —  1, ...,  N.  Hence, 

max  (|Ain|  +  \A2,n\}  <  8s2.  (10.11) 

1  <n<N 

Introduce  the  positive  constant  M*  =  M*  ^ 1 1  O'*  1 1  c2+« /n)  x [s  s]  =  M *  (C* ,  s)  by 


M* 


2(7*  max  8s2,  max  {|A1>n|  +  |A2)„|} 


l<n<N 


(10.12) 
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Hence,  (10.10)-(10.12)  imply  that 


(10.13) 


M*  <  1 6C*s2. 

Consider  the  Dirichlet  boundary  value  problem 

3 

Aw  +  bj{x)uXj  —  d(x)u  —  f  (x) ,  x  G  ft, 

3= 1 

u  |en=  9  {x)  G  C2+a  (0ft) . 

Assume  that  the  following  conditions  are  satisfied 

bj,d,f  G  Ca  (ft)  ,d(x)>  0;  max  (\bj\a,  |d|J  <  Cu 

where  C\  is  a  positive  constant.  By  the  Schauder  theorem,  there  exists  unique  solution 
u  G  C2+a  (ft)  of  this  boundary  value  problem,  and  with  a  constant  K  =  K  (Ci,ft)  >  0 
depending  only  on  the  number  C\  and  the  domain  ft  the  following  estimate  holds 

M2+C*  —  K  llfi'llc2+“(c)n)  +  \f\ot  ■ 

In  the  formulation  of  Theorem  10.1  we  provide  estimates  (10.17)-(10.22)  via  M*  and  also 
use  (10.14)  to  obtain  estimates  via  s. 

Theorem  10.1.  Let  ft  C  M3  be  a  convex  bounded  domain  with  the  boundary  d ft  G  C3. 
Suppose  that  the  inequality  (10.12)  holds.  Let  the  exact  coefficient  c*  (x)  satisfies  (2.6),  (2.7) 
and  also  c*  >  2d,  where  the  numbers  d\,d2  >  0  are  given.  For  any  function  c(x)  G  Ca  (M3) 
such  that  c  (x)  >  d\  in  ft  and  c  (x)  =  2d\  in  M3\ft  consider  the  solution  wc  (x,  s)  G 
C3  (M3\  {|x  —  xq\  <7  }),Vy  >  0  of  the  problem  (2.9a, b).  Let  Vc{x)  =  s~2  lnwc  (x,  s)  G 
C2+a  (ft)  be  the  corresponding  tail  function.  Suppose  that  the  cut-off  pseudo  frequency  s  is 
so  large  that  both  for  c*  (x)  and  any  such  function  c  (x)  the  following  estimates  hold 

\V*\2+a<f,\Vc\2+Q<f,  (10.15) 

where  (  G  (0, 1)  is  a  sufficiently  small  number.  Let  i  (x,s)  G  C2+a  (ft)  be  the  initial  tail 
function  and  let 

l^i,i|2+Q<e  (10.16) 

Denote  rj  2  (h  +  a  +  £  +  e) .  Choose  an  arbitrary  constant  C\  independent  on  other  above 
parameters.  Let  K  =  K  (C\ ,  ft)  >  0  be  the  constant  of  the  Schauder  theorem  and  N  <  N  be 
the  total  number  of  functions  qn  calculated  by  the  above  algorithm.  Suppose  that  the  number 
N  =  N  (h)  is  connected  with  the  step  size  h  via  N  ( h )  h  =  (3,  where  the  constant  fi  >  0  is 
independent  on  h.  Let  (3  be  so  small  that 


(2  1  Ci  \  .  (2  1  C7C*  \ 

y7’  16 2KC*~sA)  16 2C*~sA )  ~  \7’  16 KM*'  (M*)2 /  ’ 


(10.17) 
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In  addition,  let  the  number  77  and  the  parameter  A  of  the  CWF  satisfy  the  following  estimates 


V  <  Vo  {^,M*,d1,C1)  =  v0  [V,  lklC2+«(n)xCiM 


1  1 

=  nun  I  - 


d\ 


C 1 


2  ’  AK ’  32  •  16C*s4  ’  8C*s2  1  ~ 


^..11  di  2C'i 
<  mm  - 


2’  4 K1  32 M*s21  M*  /  ’ 


A  >  A0  ( C *,  K,  s,  v)  =  max  (  164  (C*)4  s8,  6  •  16"  (C*Y  Ks\  K 

V  n 

Then  for  every  integer  n  G  [l,lV]  the  following  estimates  hold 


?2  t r^*\ 2  jo'— 4 


|?n,fc  -  CI2+q  <  2  KM*  —  +  3V)<  32C*Ks2  -=  +  3r;  , 


Va 


v/A 


I Qn,k 1 2+a  <  2AT  <  32C*s2, 


(10.18) 

(10.19) 


(10.20) 

(10.21) 


128C*s4  (  -4=  + 

Vv/A 

/n  addition,  functions  cn ^  (. x )  >  d\  in  Q  and  cnjfc  (x)  >  d\  in  M3. 

Remarks: 

1.  It  often  happens  in  the  computational  practice  of  inverse  problems  that  theoretical 
estimates  in  convergence  theorems  are  more  pessimistic  than  ones  obtained  in  numerical 
studies.  Our  computational  experiments  show  that  this  is  exactly  our  case  in  reference  to 
estimates  (10.17)-(10.22).  For  this  reason,  we  use  in  our  computations  a  stopping  rule,  which 
is  different  from  (10.17). 

2.  By  the  Tikhonov  concept,  the  constant  C*  should  be  known  a  priori.  It  is  reasonable 
to  assume  that  C*  is  independent  on  s,  although  we  do  not  use  this  assumption  in  our  proof. 

3.  The  parameter  77  characterizes  the  error  both  in  the  data  and  in  our  mathematical 
model.  One  should  have  77  — >  0.  However,  since  in  the  reality  it  is  off  its  limiting  value  and  we 
also  have  some  other  parameters,  it  is  important  to  conduct  numerical  experiments,  which 
would  verify  this  theorem. 

4.  Truncating  integrals  at  a  high  pseudo  frequency  s  is  a  natural  thing  to  do,  because 
one  routinely  truncates  high  frequencies  in  physics  and  engineering.  By  truncating  integrals, 
we  actually  come  up  with  a  different,  although  a  quite  reasonable  mathematical  model. 
Consider  now  the  influence  of  this  truncation  on  the  accuracy  of  the  reconstruction.  Let, 
for  example  h  =  e  =  cr  =  f,  and  A-1^2  =  f.  Then  (10.22)  implies  that  the  error  of  our 
reconstruction  is  O  (£)  for  £  — »  0.  In  other  words,  one  of  claims  of  Theorem  10.1  is  that  the 
error  of  the  reconstruction  of  the  unknown  coefficient  is  mainly  determined  by  the  truncation 
error,  which  means  the  error  in  our  new  mathematical  model. 

5.  Conditions  (10.15)  and  (10.16)  with  a  small  number  £  are  natural  ones,  because  the 
number  s  is  supposed  to  be  sufficiently  large,  and  by  (3.4)  the  function  v  (x,s)  tends  to  zero 
together  with  its  ^-derivatives  as  s  — >  00.  Therefore,  the  condition  (10.16)  does  not  imply 


\Cn,k  ~  c*\a  <  8M*s  (  -j=  +  3rj  )  < 


377 


(10.22) 
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(b) 


(c) 


Figure  25:  The  hybrid  mesh  (b)  is  a  combinations  of  a  structured  mesh  (a),  where  FDM  is 
applied,  and  a  mesh  (c),  where  we  use  FEM,  with  a  thin  overlapping  of  structured  elements. 


the  assumption  of  the  closeness  of  the  first  guess  to  the  correct  solution.  Furthermore,  in 
our  computations  we  always  take  the  initial  tail  =  0,  which  reflects  the  fact  that  no  a 
priori  knowledge  about  the  medium  is  given. 

6.  One  of  the  back  bones  of  the  theory  of  ill-posed  problems  is  that  the  number  of 
iterations  can  be  chosen  as  a  regularization  parameter,  see,  e.g.,  page  157  of  [33].  Therefore, 
we  have  a  vector  (s,  N,  m\, ...,  mjfj  of  regularization  parameters,  see  details  about  their 
choice  in  subsection  7.4.  Setting  N  (h)  h  =  (3  =  const.  >  0  is  in  an  agreement  with,  e.g., 
Lemma  6.2  on  page  156  of  [33],  since  this  lemma  shows  a  connection  between  the  error  in 
the  data  and  the  number  of  iterations  (that  lemma  is  proven  for  a  different  algorithm). 

7.  It  seems  to  be  at  the  first  glance  that  because  of  (10.22),  one  can  stop  the  iterative 
process  at  n  =  1.  However,  our  numerical  experience  shows  that  this  way  one  cannot  obtain 
good  images. 

10.4  Some  numerical  results 

We  work  with  the  computationally  simulated  data.  That  is,  the  data  are  generated  by 
computing  the  forward  problem  (10.24)  with  the  given  function  c(x ).  To  solve  the  forward 
problem  (10.24),  we  use  the  so-called  hybrid  FEM/FDM  method.  The  computational  domain 
in  all  our  tests  G  =  Gfem  U  Gfdm  is  set  as  G  =  [—4.0, 4.0]  x  [—5.0,  5.0].  This  domain  is 
split  into  a  finite  element  domain  Gfem  :=  H  =  [—3.0,  3.0]  x  [—3.0,  3.0]  and  a  surrounding 
domain  Gfdm  with  a  structured  mesh,  see  Figure  25.  The  space  mesh  consists  of  triangles 
in  fl  and  of  squares  in  Gfdm-,  with  the  mesh  size  h  =  0.125  in  the  overlapping  regions.  At 
the  top  and  bottom  boundaries  of  G  we  use  first-order  absorbing  boundary  conditions,  which 
are  exact  in  this  particular  case.  At  the  lateral  boundaries,  mirror  boundary  conditions  allow 
us  to  assume  an  infinite  space  domain  in  the  lateral  direction. 
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The  forward  problem  (10.24)  is  computed  in  the  rectangle  G  C  M2  (Figure  21).  The 
coefficient  c(x)  is  unknown  only  in  a  square  Q  C  G  and 


c(x)  =  1  in  G\Q . 


(10.23) 


Hence,  2 d\  =  1  in  (2.6),  (2.7).  The  trace  of  the  solution  of  the  forward  problem  is  recorded 
at  the  boundary  dil.  This  trace  generates  the  Dirichlet  boundary  data  ij)(x,s)  in  (10.1) 
(after  the  Laplace  transform).  Next,  the  coefficient  c(x)  is  “forgotten1',  and  our  goal  is  to 
reconstruct  this  coefficient  for  x  G  Q  from  the  data  ^  (x,  s ) .  The  boundary  of  the  rectangle 
G  is  dG  =  dGiUdG2^dG$.  Here,  dG\  and  dG2  are  respectively  top  and  bottom  sides  of  the 
largest  rectangle  of  Figure  25  and  dG 3  is  the  union  of  left  and  right  sides  of  this  rectangle. 
The  forward  problem  is 

0,  in  G  x  (0,  T), 

Qu 

°,  ^r(-,°)  =  0,  in  G, 

f  (t) ,  on  dGi  x  (0,  H],  (10.24) 

dtu,  on  dG\  x  (H,  T), 
dtu,  on  dG2  x  (0,  T), 

0,  on  d G3  x  (0,  T), 

where  T  is  the  final  time.  When  calculating  the  Laplace  transform  of  the  boundary  data, 
we  integrate  for  t  G  (0,  T),  thus  calculating  an  approximation  of  this  transform.  The  plane 
wave  /  is  initialized  at  the  top  boundary  dGi  of  the  computational  domain  G,  propagates 
during  the  time  period  (0,H]  into  G ,  is  absorbed  at  the  bottom  boundary  dG2  for  all  times 
t  G  (0,  T)  and  it  is  also  absorbed  at  the  top  boundary  dGi  for  times  t  G  (G,  T ).  Here 


,  .  o  u 

c  (x)  ^777  —  A u  = 


dt 2 


u(-,0)  = 

^nU\dG1  = 
^nU\dGi  = 
®nU\dGi  = 
®nU\  dG3  = 


m 


(sin  (st— 7t/2)  +  1) 


10 


5 


0  <  t  <  H  := 


17.8H. 


10.4.1  Main  discrepancies  between  the  theory  and  the  numerical  implementa¬ 
tion 

It  makes  sense  to  summarize  in  this  subsection  main  discrepancies  between  the  above  theory 
and  our  numerical  implementation.  We  note  that  such  discrepancies  quite  often  occur  in 
computations  of  inverse  problems. 

The  first  discrepancy  is  that  in  order  to  generate  the  data  for  the  inverse  problem,  we 
have  solved  the  forward  problem  (10.24)  in  the  finite  domain  G  with  the  plane  wave  instead 
of  the  Cauchy  problem  (2.1),  (2.2)  with  the  point  source.  To  explain  this,  we  note  that 
our  theory  needs  the  problem  (2.1),  (2.2)  only  for  the  asymptotic  behavior  (2.11),  and 
consequently  for  (2.16).  However,  conditions  guaranteeing  the  asymptotic  behavior  (2.11) 
(see  those  conditions  in  [12])  are  very  hard  to  verify  computationally.  Therefore,  in  all  our 
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numerical  tests  we  have  somewhat  verified  numerically  the  asymptotic  behavior  (2.16).  To 
do  this,  we  have  considered  functions  g\  (s)  and  r/2  (s) ,  for  s  E  [6.5,  7.5]  D  [s,  s]  =  [6.7,  7.45] , 
where 

9i  (s)  =  s  ||  Vv  (x,  s )  || ia(n) ,  g2  (s)  =  s2 1|  Vg  {x,  s )  || La(n) , 

where  functions  v  (x,  s)  and  q  (. x ,  s)  are  taken  from  the  solution  of  the  forward  problem.  We 
have  found  that  the  interval  [s,  s]  =  [6.7,7.45]  was  the  optimal  one  for  our  reconstruction 
procedure,  and  we  have  used  it  in  all  our  numerical  experiments.  Graphs  of  functions  g±  (s) 
and  r/2  (s)  (not  presented  here)  have  shown  that  these  functions  are  very  close  to  constants 
for  s  E  [6.5,  7.5] ,  which  corresponds  well  with  (2.16).  We  have  also  verified  numerically  in 
all  our  tests  that  the  function  w(x,s)  >  0,  which  justifies  the  introduction  of  the  function 
v(x,s)  =  lnw(x,s)  / s 2.  In  our  opinion,  these  verifications  provide  numerical  justifications 
for  the  above  discrepancy. 

We  describe  the  2  nd  main  discrepancy  in  this  paragraph.  Because  of  (10.23),  we  set 
2di  =  1.  Instead  of  using  the  extension  procedure  described  in  the  beginning  of  section  5, 
we  simply  set  cUtk  (x)  :=  1  in  G\lf.  In  addition,  since  by  (2.6)  we  need  a  priori  lower  bound 
c(x)  >  di,  we  enforce  that  the  coefficient  c(x)  belongs  to  the  set  of  admissible  coefficient 
Cadm  =  {c(x)  >  0.5}  as  follows:  If  cn^(x o)  <  0.5  for  a  certain  point  xq  E  If  and  a  certain  pair 
(■ n,k ),  then  we  set  cn^(x o)  :=  1.  The  only  reason  why  we  use  the  value  1  in  this  setting  is 
that  we  are  supposed  to  know  that  the  condition  (10.23)  is  satisfied.  Therefore,  this  setting 
as  well  as  the  fact  that  we  allow  the  function  c(x)  to  attain  values  between  0.5  and  1  does 
not  mean  that  we  assume  the  knowledge  of  the  background  value  of  the  function  c(x).  In 
principle,  this  constraint  cannot  guarantee  neither  the  continuity  of  the  resulting  function 
cn,k  (x)  nor  that  cny.  (x)  >  1.  Nevertheless,  we  have  observed  in  our  numerical  tests  that 
all  resulting  functions  cn^  are  continuous  and  cn^  (x)  >  1  for  all  x,  i.e.,  “allowed”  values 
between  0.5  and  1,  are  not  actually  attained  in  iterations.  Another  reasonable  setting  would 
be  to  assign  cn^(x 0)  :=  0.5.  We  have  not  done  this  yet  and  hope  to  investigate  this  topic  in 
the  future. 

The  3rd  main  discrepancy  is  that  our  square  does  not  have  a  smooth  boundary,  as  it  is 
required  by  the  Schauder  theorem.  Furthermore  the  Dirichlct  boundary  value  problems  for 
functions  qn^  in  the  square  If  were  solved  by  the  FEM,  in  which  the  same  finite  elements  were 
used  as  ones  in  the  forward  problem  (10.24).  The  “inverse  crime”  was  not  committed  because 
the  forward  problem  was  solved  for  the  hyperbolic  equation  (10.24),  whereas  we  solve  many 
different  elliptic  equations  in  our  iterative  algorithm.  The  FEM  cannot  guarantee  that 
resulting  functions  qnx,  E  C2+a  (if)  ,  as  it  is  required  by  Theorem  10.1.  Nevertheless,  an 
analogue  of  Theorem  10.1  can  be  proved  for  the  discrete  case  when  the  FEM  analogues 
of  equations  for  functions  qntk  are  used,  and  also  the  domain  If  with  (91f  E  C3  is  replaced 
respectively  with  either  a  rectangular  prism  in  M3  or  a  rectangle  in  M2,  as  in  our  numerical 
examples.  To  prove  this  analogue,  one  needs  to  use  the  weak  formulations  of  equations  (5.4), 
(5.6)  and  the  Lax-Milgram  theorem  instead  of  the  Schauder  theorem.  Next,  because  of  the 
equivalency  of  norms  in  finite  dimensional  spaces,  the  rest  of  the  proof  is  very  similar  with 
the  above.  However,  the  latter  development  is  outside  of  the  scope  of  this  publication  and 
might  be  considered  in  our  future  works.  Another  interesting  question  here  is  about  the 
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change  of  reconstructed  images  due  to  the  increase  of  the  number  of  finite  elements,  because 
that  equivalency  of  norms  “worsens”  with  the  increase  of  the  dimension  of  the  space.  This 
question  might  also  be  addressed  in  future  publications. 

The  4th  main  discrepancy  is  that  because  of  the  above  mentioned  equivalency  of  norms, 
as  well  as  because  the  L2  (hi)  norm  is  easier  to  verify  computationally  than  the  Ca  (hi)  norm 
(especially  because  a  G  (0,1)),  we  use  L2  (hi)  norms  for  our  stopping  rule.  An  additional 
discrepancy  is  that  in  order  to  produce  updates  for  tails,  we  have  solved  on  each  iterative 
step  the  forward  problem  (10.24)  instead  of  solving  the  elliptic  problem  for  the  function  w. 
Next,  we  have  calculated  the  Laplace  transform  to  obtain  the  function  wn^  (x,s) . 


10.4.2  Results  of  reconstruction 

In  this  subsection  we  present  results  of  our  reconstructions.  We  have  performed  numerical 
experiments  to  reconstruct  the  medium,  which  is  homogeneous  with  c  (. x )  =  1  except  of 
either  two  small  squares  or  a  single  square,  see  Figure  25.  However,  we  have  not  assumed 
a  priori  knowledge  of  neither  the  structure  of  this  medium  nor  of  the  background  constant 
c(x)  —  1  outside  of  those  squares. 

In  all  our  numerical  experiments  we  have  chosen  the  step  size  with  respect  to  the  pseudo 
frequency  h  =  0.05.  Hence,  N  =  15  in  our  case.  It  is  important  that  in  all  our  tests  the 
regularization  parameters  where  chosen  the  same.  This  means  that  results  were  not  “tilted” 
towards  desired  ones.  We  have  chosen  two  sequences  of  regularization  parameters  A  :=  An 
and  £  =  en  for  n  =  1,  ...,IV.  The  formulation  of  Theorem  10.1  remains  almost  unchanged 
for  this  case.  The  reason  of  choosing  different  values  of  An  and  en  is  that  values  gradients 
of  functions  q\  and  q2  are  very  small.  Hence,  in  order  not  to  eliminate  totally  the  influence 
of  the  nonlinear  term  (V qn^~ if  ,  n  —  1,2  in  (10.8),  the  values  of  Ai  and  A2  should  not  be 
too  large.  Next,  the  values  of  the  nonlinear  term  start  to  grow,  and  we  balance  them  by 
taking  a  larger  value  of  An  for  n  =  3, 4,  5.  For  n  >  5  the  values  of  the  nonlinear  term  become 
even  bigger,  and  we  balance  them  via  increasing  the  value  of  An  again.  This  points  towards 
the  importance  of  the  introduction  of  CWFs  in  the  numerical  scheme,  as  compared  with  the 
decrease  of  the  step  size  h.  The  considerations  for  choosing  different  values  of  en  are  similar. 
In  all  our  tests  of  []  the  values  of  the  parameters  An  and  en  were: 

Xn  =  20,  n  =  1,  2;  An  =  200,  n  =  3, 4,  5;  An  =  2000,  n  >  6; 
en  —  0,  n  =  1,  2;  £n  =  0.001,  n  =  3, 4,  5;  £n  =  0.01,  n  =  6,  7, 

£n  =  0.1,  n  >  8. 


We  use  the  following  approximation  to  find  c(x)  at  the  point 


i  7  G+ 1  j  2iy?  T  ij  Vij- )-i  'Xvj  j  T  V{  j—  \ 

C  —  - — - 1 - 


dx2 


dy 2 


dx 
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-'ij+1  Ui,j 
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where  dx  and  dy  are  grid  step  sizes  of  the  discrete  finite  difference  mesh  in  the  directions 
x  and  y  respectively.  We  also  use  the  smoothness  indicator  in  values  of  cn,k(x )  via  local 
averaging  over  the  neighboring  elements. 

The  resulting  computed  function  is  c  (x)  :=  Cj y(x).  Recall  that  the  number  of  iterations  is 
a  part  of  the  vectorial  regularization  parameter  in  our  case.  We  have  used  mn  =  4  iterations 
with  respect  to  tails  for  n  <  no  and  mn  =  7  for  n  —  n0  +  1, ...,  TV,  where  numbers  no  and 
N  are  chosen  on  the  basis  of  an  objective  stopping  rule  described  below.  Hence,  while  the 
pairs  ( no,N )  differ  in  our  tests,  the  rule  of  their  choice  (i.e.,  the  stopping  rule)  remains 
the  same.  As  it  is  always  the  case  in  ill-posed  problems,  the  choice  of  proper  regularization 
parameters  and  of  a  proper  stopping  rule  was  time  consuming.  However,  we  point  out  that 
our  stopping  rule  as  well  as  regularization  parameters  Xn,  en,  mn,  s  once  chosen,  remained  the 
same  for  all  our  numerical  experiments  described  in  Tests  1-4  below.  Hence,  results  were 
not  “conveniently  adjusted”  for  each  specific  test  in  order  to  obtain  the  best  possible  image 
for  that  test. 

In  all  our  tests  we  have  introduced  the  multiplicative  random  noise  in  the  boundary  data. 
Next,  we  apply  the  Laplace  transform  to  the  boundary  data,  which  helps  to  both  “smooth 
out”  and  decrease  the  noise,  due  to  the  integration.  Because  of  that,  we  have  successfully 
used  the  following  formula  for  the  s— derivative  of  the  boundary  data  <p  (x,s)  in  (10.1) 

dip  (x,  sn)  _  ip  (x,  s„_ 1)  -  ip  (x,  sn)  ,  _  n  nK 
Os  ~  h  ,  h  —  U.U5. 

In  all  tests  of  [12]  the  starting  value  for  the  tail  was  chosen  (x,s)  =  0,  which  is  in 
an  agreement  with  (10.16)  and  also  reflects  the  fact  that  no  advanced  knowledge  about  tails 
is  available.  While  we  present  results  of  only  two  tests  here,  more  numerical  results  can 
be  found  in  [].  It  is  also  important  that  in  all  our  tests  we  have  used  the  same  objective 
stopping  rule,  which  was  based  on  stabilization  of  relative  L2  norms  of  discrepancies  between 
two  consecutive  iterations  of  both  tails  and  functions  cny.,  see  details  in  [12]. 

Test  1.  We  now  test  our  numerical  method  on  the  reconstruction  of  the  structure  given 
on  Figure  25c.  We  introduce  0  =  5%  of  the  multiplicative  random  noise  in  the  boundary 
data.  We  take  c  =  4  for  both  small  squares  of  Figure  25c  and  c  =  1  outside  of  these  squares. 
Hence,  the  inclusion/background  contrast  is  4  :  1.  Figure  26  presents  resulting  image  of  the 
function  c  :=  012,7.  Figure  27  displays  the  one-dimensional  cross-sections  of  the  images  of 
the  function  cn,fc  along  the  vertical  line  passing  through  the  middle  of  the  left  small  square. 
The  imaged  function  c(x)  is  superimposed  with  the  correct  one.  One  can  see  that  the  value 
of  the  function  c  (x)  both  inside  and  outside  of  the  inclusion  is  imaged  correctly,  although 
the  location  of  the  inclusion  is  somewhat  shifted  to  the  top. 

Test  2.  While  our  method  does  not  require  a  good  a  priori  guess  about  the  solution, 
the  main  point  of  this  test  is  to  show  that  a  reconstruction  algorithm,  which  is  based  on  the 
minimization  of  a  least  squares  objective  functional,  might  lead  to  a  poor  reconstruction,  if 
a  good  first  guess  about  the  solution  is  unavailable.  We  use  the  reconstruction  algorithm 
described  in  [26-28] ,  where  the  inverse  problem  is  formulated  as  an  optimal  control  problem 
of  the  minimization  of  a  least  squares  objective  functional.  The  latter  is  solved  by  the 
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Figure  26:  Computed  image  for  Test  1. 


Figure  27:  Computed  image  for  Test  2 


(a) 


(b) 


(c) 


Figure  28:  Computed  image  for  Test  5. 


quasi-Newton  method,  which  is  known  to  be  a  good  method  for  this  purpose.  We  find  a 
stationary  point  of  a  Lagrangian,  using  the  forward  wave  equation  (the  state  equation),  the 
backward  wave  equation  (the  adjoint  equation),  and  an  equation,  which  reflects  the  fact  that 
the  gradient  with  respect  to  the  parameters  should  vanish.  A  minimizer  of  a  corresponding 
least  squares  objective  functional  is  found  via  an  iterative  procedure  via  solving  the  forward 
and  backward  wave  equations  for  each  iterative  step  and  updating  the  unknown  coefficients. 
We  generate  the  data  for  the  inverse  problem  using  the  same  computational  mesh  and  the 
same  parameters  as  ones  in  Test  1.  We  start  the  optimization  algorithm  with  different  values 
of  the  first  guess  for  the  parameter  cguess  =  const,  at  all  points  of  the  computational  domain 
fh  Figure  28  presents  the  images  of  the  computed  function  ccomp.  for  the  following  initial 
guesses:  on  a)  cguess  =  1.0,  on  b)  cguess  =  1.5,  and  on  c)  cguess  =  2.0.  We  observe  that  images 
deteriorate  from  a)  to  c)  with  the  deterioration  of  the  first  guess.  We  conjecture  that  local 
minima  are  achieved  in  all  these  three  cases,  although  an  investigation  of  the  latter  topic  is 
outside  of  the  scope  of  this  publication. 

We  note  that  in  [26-28]  the  quasi-Newton  method  was  applied  in  a  combination  with 
the  so-called  “adaptivity  technique” .  ft  is  known  that  this  technique  is  capable  to  improve 
the  quality  of  some  images,  which  are  initially  obtained  via  minimizations  of  least  squares 
objective  functionals.  We  are  not  presenting  corresponding  results  here  because  of  space 
considerations.  We  refer  to  the  recent  work  [28]  where  tests  were  conducted  using  the 
adaptivity  technique  for  a  structure  similar  with  one  considered  in  Test  2.  It  was  concluded 
in  [28]  that  while  the  adaptivity  significantly  improves  images  with  the  cguess  =  1.0  and 
cguess  =  1-5,  the  image  quality  with  cguess  =  2  has  deteriorated  not  only  on  the  coarse  mesh 
but  also  on  the  one  and  two  times  refined  meshes.  Hence,  the  adaptivity  technique  works 
well  in  [28]  only  in  a  neighborhood  of  the  initial  guess  cguess  G  [1, 1.5]  . 
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10.4.3  An  extension  of  the  method  of  [12]. 

In  the  joint  effort  of  the  PI  with  Drs.  H.  Shan,  J.Su,  N.  Pantong  and  H.Liu,  who  were 
co-authoring  in  [8],  an  extension  of  the  above  idea  of  this  section  for  the  case  of  the  running 
source  was  carried  out  [8].  Since  there  is  no  clear  asymptotic  behavior  of  the  solution  of  the 
equation  ()  when  the  source  runs  to  infinity,  a  heuristic  treatment  of  tails  was  proposed  in 
[8].  Images  of  [8]  are  of  the  same  quality  as  ones  of  section  6,  see  Figures  15  and  16.  In 
addition,  recently  we  have  figured  out  how  to  arrange  a  “clean”  asymptotic  behavior  and  this 
work  is  currently  in  progress.  Preliminary  results  on  the  latter  idea  are  encouraging. 


11  Some  Analytical  Results  of  the  Project 


Although  the  main  target  of  this  project  is  globally  convergent  numerical  methods  for  CIPs, 
it  is  quite  natural  that  a  number  of  purely  analytical  results  were  also  obtained,  in  addition 
to  the  above  uniqueness  Theorem  8.1.  These  results  are  reflected  in  publications  [15-20].  In 
this  section  we  outline  most  interesting  analytical  results,  also  see  annual  report  of  2006  for 
more  details  [21]. 

11.1  Inverse  problems  for  the  non-stationary  transport  equation 

Results  of  this  subsection  have  formed  a  major  part  of  the  Ph.D.  thesis  of  Mr.  Sergey  E. 
Pamyatnykh  [15],  [16],  whose  Thesis  Advisor  was  the  PI.  Sergey  has  defended  his  thesis  in 
2006  at  The  University  of  North  Carolina  at  Charlotte.  No  analogs  of  Theorem  11.1  were 
known  before  the  publication  [15].  Theorem  11.2  [16]  is  the  first  global  uniqueness  result  for 
a  non-overdetermined  CIP  for  the  transport  equation.  Previous  uniqueness  results  for  CIPs 
for  the  transport  equation  were  obtained  either  for  the  case  of  the  over-determined  data 
or  under  some  restrictions  imposed  on  coefficients  of  the  transport  equation,  see  [23,48,49]. 
This  is  because,  unlike  results  of  this  subsection,  the  method  of  Carleman  estimates  was  not 
applied  previously. 

Let  R  =  const.  >  0, 

n  =  {x  E  :  |z|  <  R},  Sn  =  {ueRn  :  \u\  =  l},Z  =  Ox  Sn , 

H  =  Q,  x  (— T,  T)  xSn,  r  =  <9f lx  (— T,  T)  x  Sn, 

H+  =  n  x  (0,  T)  x  Sn,  r+  =  dQ  x  (0,  T)  x  Sn, 

H~  =  flx  (— T,  0)  x  Sn,  r~  =  dn  x  (— T,  0)  x  sn, 

so  that 


h  =  h+uh-  and  r  =  r+ur~. 
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Also,  introduce  spaces 


Ck(H),  k  =  1,2 ,Ck(H)  =  | u(x,t,is)  :  maxJ\Dsxtu(x,t,is)\\c(S)  =  \\u\\dk(n)  <  ooj  . 

The  non-stationary  transport  equation  in  the  domain  H  has  the  form  [6] 

Ku  :=  ut  +  (is,  Vm)  +  a(x,  t,  is)u  H —  J g(x,  t,  is,  p)u(x,  t,  p)du^  =  F(x,  t,  is),  (11.1) 

sn 

where  v  G  Sn  is  a  unit  vector  of  the  particle  velocity,  u(x,t,v)  G  C1(i/)  is  the  density  of 
the  particle  flow,  a(x,t,v)  is  the  attenuation  coefficient,  F(x,t,is )  is  the  angular  density  of 
sources,  g(x,t,is,n )  is  a  scattering  indicatrix.,  and  (i/,  Vm)  denotes  the  scalar  product  of 
vectors  is  and  Vw  in  Mn. 

Consider  the  following  boundary  condition 

u(x,  t ,  is)  =  p(x,  t ,  is),  (11.2) 

for  {x,  f,i/)ern  {(n(x),  is)  <  0}  . 

Here  (n(x),  is)  is  the  scalar  product  of  the  outer  unit  normal  vector  n{x)  on  the  surface  dVt 
and  the  direction  of  the  velocity  is.  Hence,  only  incoming  radiation  is  given  at  the  boundary 
in  this  case.  The  boundary  condition  (11.2)  together  with  the  initial  condition  at  t  =  —T 

u(x,—T,is)  —  f(x,v),  V(x,v)  £  Z  (11.3) 

form  a  classic  boundary  value  problem  for  the  non-stationary  transport  equation  (11.1)  in  the 
domain  H.  This  problem  means  that  given  the  initial  condition  and  the  incoming  radiation 
at  the  boundary,  find  the  density  of  the  particle  flow. 

Consider  now  a 

A  non-standard  boundary  value  problem.  Suppose  that  the  initial  condition  (11.3) 
is  unknown,  but  the  outgoing  radiation  at  the  boundary  is  known  instead.  In  other  words, 
the  following  function  h(x,t ,  is)  is  known  in  addition  to  the  function  p(x,t,  is)  in  (11.2), 

u(x,t,is)  —  h(x,t,is),  (11.4) 


for  (x,t,v)  G  T  fl  {(n(x),  is)  >  0}  . 

Determine  the  function  u(x,t ,  is)  in  the  domain  H. 

The  problem  (11.1),  (11.2),  (11.4)  is  ill-posed.  Hence,  one  cannot  expect  to  obtain  an 
existence  theorem.  However,  one  can  expect  to  obtain  a  stability  result,  i.e.,  to  estimate  the 
function  u  via  functions  p  and  h. 
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Theorem  11.1  (Lipschitz  stability)  [15].  Let  in  (7.1)  ||a||c^  <  r  and  ||g||cp 


(. HxSn )  —  r’ 

where  r  =  const.  >  0.  Assume  also  that  the  function  F  e  L2  (H)  and  T  >  R.  Denote 


Q  (x,  t,  v)  = 


p(x,t,is)  if  (n(x),is)  <  0 
h(x,t,v )  if  (n(x),is)  >  0 


Let  the  function  u  G  C1  (12  x  [— T,  T])  x  C  (S'")  satisfies  conditions  (6. If),  (6.15),  (6.17). 
Then  the  following  Lipschitz  stability  estimate  is  valid 


Ml l2(h)  ^  C 


+  Mil 


L2(H) 


where  the  positive  constant  C  =  C  (■ r ,  R,  T )  depends  only  on  numbers  r,  R  and  T . 

The  PI  and  Pamyatnykh  [16]  have  considered  the  following 

Inverse  Problem.  Suppose  that  the  equation  (11.1)  is  satisfied  in  the  domain  H+  = 
12  x  (0,  T )  x  Sn.  Suppose  that  the  attenuation  coefficient  a  (. x ,  is)  is  independent  on  time  t 
and  unknown,  but  the  following  functions  s  ( x ,  is)  and  k  (x,  t,  is)  are  known 


(11.5) 

(11.6) 


C(H+xSn)  — 


U  (x,  0,  is)  —  S  (x,  is) , 
u  (. x ,  t,  is)  |r+=  k  (x,  t,  is) . 

Determine  the  function  a  (x,  is) . 

Thus,  the  following  theorem  was  proven 

Theorem  11.2  [16].  Suppose  that  the  derivative  dtg  exists  in  H+xSn  and 
r i  for  k  =  0,1,  where  rq  is  a  positive  constant.  Let  in  (7.6)  |s(x,  i>)|  >  r2,  where  r2  = 
const  >  0.  Assume  that  there  exist  two  pairs  of  functions  ( a\,u\ )  and  (a2,u2)  satisfying 
(11.1),  (11.2),  (11.5),  (11.6)  and  such  that 

«1,  a2  G  C  ( Z )  and  uu  uit,  uitt,  Vuu  Vuit  G  C  (. H+ )  ,  i  =  1,  2. 

Suppose  also  that 

[(dif)  (x,  is)}2  =  [(a,;/)  (x,  -is)]2  ,i  =  1,2. 

Let  <  r3,  where  r3  =  const  >  0.  Then  there  exists  a  number  T0  =  T0  (R,  r2,r3)  > 

R  such  that  if  T  >  T0,  then  a±  —  a2  in  Z  and  v.\  =  m2  in  H+ .  If  the  function  m2(x,  t,  is)  ^  0 

in  H+,  then  it  is  sufficient  to  assume  that  T  >  R. 

11.2  Exact  controllability  for  the  non-stationary  transport  equa¬ 
tion 

In  the  work  of  Klibanov  and  Yamamoto  [17]  the  exact  controllability  theorem  for  the  non- 
stationary  transport  equation  was  proved  for  the  first  time.  Consider  the  homogeneous 
transport  equation  (11.1), 


Ku  :=  ut  +  (is,  Vu)  +  a(x,  t ,  is)u  +  /  g(x,  t,  v,  n)u(x,  t ,  pTjdcr^  =  0  in  H+, 


(11.7) 
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where  H+  is  defined  above  with  the  only  difference  that  now  12  C  Mn  is  a  bounded  domain 
with  its  boundary  <912  G  C°° .  Assume  that  functions 


aeC1  (H+)  ,geCl(H+  x  Sn )  . 


(11.8) 


Let 

=  {(x,  t,  v)  G  <912  x  (0,  T)  x  Sn  :  (n(x),  v)  <  0}  , 

Introduce  the  Hilbert  space  of  real  valued  functions  L2cos  (r_)  as  the  one  with  the  scalar 
product 

(p,  q)  =  j p  (x,t,u)  q  (x,t,u)  |cos  (n(x),  v)\  dSxdtdau. 
r_ 

Observe  that  a  particular  difficulty  in  working  with  this  scalar  product  is  that  the  weight 
function  |cos  ( n(x ),  v)\  has  zeros  on  r_. 

We  first  introduced  the  weak  solution  u  E  C  ([0,  T]  (Ox  Sn ))  of  the  equation  (11.7) 
with  the  initial  condition 

u  (x,  0,  v)  —  0  (11.9) 

and  the  boundary  condition 

u  |r„=  P  (x,  t,  is)  G  L2os  (r_) .  (11.10) 

To  do  this,  we  have  used  modified  classic  density  arguments.  Next,  we  posed  the  following 
The  Exact  Controllability  Problem.  Consider  an  arbitrary  function  w  (x,  v)  G 
L 2  (O  x  Sn ) .  Find  such  a  control  function  pw  (x,  t,  v)  G  L2cos  (T_)  that 

u  (x,  T,  v)  —  w  (x,  u) ,  (11.12) 

where  u  G  C  ([0,  T] ;  L2  (H  x  Sn))  is  the  weak  solution  of  the  equation  (11.7)  with  the  initial 
condition  (11.9)  and  the  boundary  condition  (11.10),  in  which  p  =  pw .  The  time  T  is  the 
so-called  “steering  time” . 

We  have  modified  Theorem  11.1  to  prove 

Theorem  11.3  [17].  Let  12  be  strictly  convex  bounded  domain  with  <912  G  C°°  and 
conditions  (11.8)  hold.  Let  the  steering  time  T  >  diameter  (12) .  Then  for  any  function 
w  (x,  v)  G  L2  (12  x  Sn )  there  exists  a  control  function  pw  (x,  t,  v)  G  L2os  (T_)  such  that  if  the 
function  u  G  C  ([0,  T] ;  L2  (12  x  Sn ))  is  the  weak  solution  of  the  equation  (11.7)  satisfying 
(11.9)  and  (11.10)  with  p  =  pw  in  (11.10),  then  (11.12)  holds. 

11.3  Estimates  of  initial  conditions  of  parabolic  equations  and  in¬ 
equalities  with  the  lateral  Cauchy  data  in  finite  domains 

The  result  of  this  subsection  is  published  in  the  work  of  Klibanov  [18].  As  it  was  stated  in 
subsection  1.2,  this  results  was  featured  as  one  of  the  best  of  2006  by  the  Editorial  Board  of 
Inverse  Problems.  Let  12  C  be  a  bounded  domain  with  the  boundary  <912  G  C1.  For  any 
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T  =  const.  >  0  denote  Qt  =  fl  x  (0,  T) ,  St  =  dfl  x  (0,  T) .  For  any  function  s(x),x  G  Mn 
denote  Si  =  ds/dxi,i  =  1, n,  whenever  the  differentiation  is  appropriate.  We  also  denote 
Vs  =  (si, sn) .  Let  L  =  L(x,  t,  D )  be  an  elliptic  operator  of  the  second  order  in  QT, 

n  n 

Lu  :=  L(x,  f,  D)u  =  a13  (x,  t)uij  +  IP  (x,  t)uj  +  b°(x,  t)u, 

i,j= i  i,j= i 


with  its  principal  part  L0, 


Lqu  :=  L0(x,t,  D)u  =  a23  (x,t)uij, 

i,j= 1 


where  coefficients 


ai3  =  a31,  aij  G  C1  (Qt)  n  5  (Qt)  ;  6°  G  B  (Qr)  , 

where  i?  (Qr)  is  the  set  of  functions  bounded  in  QT.  Naturally,  we  assume  the  existence  of 
two  positive  numbers  dq,  <72,  04  <  dq  such  that 

n 

l?|2  <  E  <  ^2  l^l2 ,  V(x,t,£)  G  Qr  x  ML  (11.13) 

*3=1 

Let  the  function  u  G  FT2,1  (Qt)  be  a  solution  of  the  parabolic  equation 

Ut  =  Lu  +  / (x,  t),  a.e.  in  Qt,  (11.14) 

with  the  unknown  initial  condition  g(x), 

u(x,  0)  =  g(x)  G  H 1  (hi) ,  (11.15) 

where  the  function  /  G  L2  (Qt)  ■  Along  with  the  equation  (11.14)  we  also  consider  a  more 
general  case  of  the  parabolic  inequality 

/  (lit  ~  Lu)2  dxdt  <  M2,  (11.16) 

Qt 

where  the  function  u  G  H 21  (Qt)  satisfies  condition  (11.15)  and  M  =  const.  >  0. 

Inverse  Problem.  Assume  that  the  following  lateral  Cauchy  data  h\  (x.  t)  and  h^ix. t) 
are  given 

du 

u\St=  h(x,t),  —  |St=  h2(x,t),ST  =  S  x  (0,T) ,  (11-17) 

where  the  function  u  G  H 2,1  (Qt)  satisfies  either  conditions  (11.14)  and  (11.15)  or  condi¬ 
tions  (11.15)  and  (11.16).  Estimate  the  unknown  initial  condition  g  and  the  function  u  via 
functions  hi,h,2  and  /. 
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This  is  an  inverse  problem  of  the  determination  of  the  initial  condition  in  the  parabolic 
equation  using  lateral  Cauchy  data  hi  and  h 2.  Applications  are  in  such  diffusion  and  heat 
conduction  processes,  in  which  one  is  required  to  determine  the  initial  state  using  boundary 
time  dependent  measurements.  To  describe  a  more  specific  application,  consider  the  cooling 
process  of  a  solid,  which  is  contained  in  a  bounded  domain  fl  C  M3.  Suppose  that  the  initial 
temperature  of  this  solid  is  high,  unknown,  and  is  the  subject  of  one’s  interest.  Suppose  also 
that  interior  points  of  this  solid  are  unavailable  for  the  temperature  measurements.  Instead 
one  is  measuring  the  time  dependence  of  both  the  temperature  u  and  the  heat  flux  at  the 
boundary  <9fi.  Assuming  that  near  <9fl  the  principal  part  L0  =  A,  we  obtain  that  the  heat 
flux  at  the  boundary  is  du/dn  |gT  .  Hence,  in  this  application  the  Inverse  Problem  is  the 
problem  of  the  determination  of  the  spatial  distribution  of  the  initial  temperature  u(a:,0), 
using  time  dependent  boundary  measurements. 

A  particular  benefit  of  considering  this  applied  example  is  that  it  helps  to  understand 
the  naturality  of  imposing  a  priori  bound  on  the  L2  (fl)  —norm  of  the  gradient  || 
in  Theorem  11.4.  Indeed,  this  assumption  means  a  priori  knowledge  of  the  absence  of  high 
gradients  in  the  initial  temperature,  which  is  quite  natural  in  this  application.  A  similar 
idea,  although  in  a  more  general  form,  is  one  of  building  blocks  of  the  theory  of  ill-posed 
problems  and  it  follows  from  the  above  mentioned  fundamental  Tikhonov  theorem  [50]. 

The  idea  of  the  proof  of  Theorem  11.4  (below)  is  to  combine  two  types  of  Carleman 
estimates:  “lateral”  and  “backwards”  ones.  Lateral  Carleman  estimate  is  the  one,  which 
estimates  the  solution  of  the  parabolic  equation  (or  inequality)  via  lateral  Cauchy  data. 
It  estimates  that  solution  in  a  subdomain  G  of  the  time  cylinder  Qt ■  However,  since  G  fl 
{t  =  0}  =  0,  then  the  lateral  Carleman  estimate  does  not  provide  an  estimate  for  the  initial 
condition  g(x).  Still,  it  ensures  an  estimate  of  the  norm 

IMMo)||L2(q)  (11.18) 

via  norms  of  the  lateral  Cauchy  data.  In  (11.18)  t0  e  (0,  T)  is  a  certain  constant. 

Thus,  one  should  somehow  estimate  the  initial  condition  g(x )  via  the  norm  (11.18).  To 
do  so,  a  backwards  Carleman  estimate  is  used.  The  backwards  Carleman  estimate  is  the  one, 
which  estimates  the  solution  u(x,t )  of  the  parabolic  equation  for  t  e  (0,  to)  0  (0,  T)  via  the 
norm  (11.18),  i.e. ,  it  estimates  solutions  of  parabolic  equations  with  the  reversed  direction 
of  time.  However,  the  previously  known  such  estimate  enabled  one  to  estimate  the  certain 
norms  of  the  function  u(x,  t )  only  in  H  x  (r,  t0 )  for  a  0  <  r  <  t0.  Hence,  the  case  of  u(x,  0)  is 
a  more  delicate  one.  The  main  new  element  of  the  work  [18]  is  a  new  backwards  Carleman 
estimate,  which  enabled  to  estimate  \\u(x,  0)||i;^^ ,  as  well  as  to  consider  the  general  case  of 
the  operator  L  with  (x,  t)—  dependent  coefficients,  including  the  inequality  (11.16). 

A  similar  problem  was  considered  earlier  by  Xu  and  Yamamoto  [51].  However,  their 
technique  cannot  handle  the  case  of  non-self  adjoint  operator  L  with  ( x ,  t )  —  dependent 
coefficient,  including  the  inequality  (11.16).  This  is  because  they  use  the  so-called  “loga¬ 
rithmic  stability”  method  for  the  backwards  estimate,  and  this  method  works  only  with  the 
self-adjoint  operator  L  with  rc-dependent  coefficients,  see,  e.g.,  [22]  for  that  method. 

The  following  stability  estimate  was  proven  in  [18]. 
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Theorem  11.4.  Assume  that  the  above  conditions  imposed  on  the  coefficients  of  the 
elliptic  operator  are  fulfilled.  In  the  case  of  the  problem  (11.15)-(11.1 7)  denote  f(x,t)  =  M. 
In  both  problems  (11. If),  (11.16),  (11.17)  and  (11.15)-(11.17)  denote  F=  (hi,h2,f)  and 

11-^11  =  Ill'll  I//1  (St)  +  IIMI  l2(st)  +  II/IIl2(Qt) 

Suppose  that  ||F||  <  B,  is  a  positive  number.  Assume  that  the  function  g  G  H1  (12) .  Then 
there  exists  a  positive  constant  C  such  that,  for  every  number  a  G  (0,  2)  there  exists  a  number 
(3  G  (0, 1)  such  that  the  following  stability  estimate  holds 

Nll,n,  <  l  IIVsIIU  ■  1"  [jjpjf]  1+c(f)“l|f’ir-“.  (11-19) 

The  constant  C  depends  only  on  the  operator  L  and  the  domain  Qt  and  the  constant  (3 
depends  on  the  same  parameters  as  C,  as  well  as  on  a. 

The  estimate  (11.19)  is  the  so-called  conditional  stability  estimate,  because  it  assumes 
boundedness  of  a  stronger  norm  ||  V^||^2^,  see  above  about  the  applied  aspect.  Conditional 
stability  estimates  are  typical  ones  for  ill-posed  problems,  as  it  follows  from  the  Tikhonov 
theorem  [50].  While  constants  like  B,  a  and  (3  do  not  appear  in  traditional  well-posed 
problems  (at  least  in  the  linear  case),  their  appearance  is  quite  natural  in  ill-posed  problems. 
On  the  other  hand,  assuming  that  ||F||  is  sufficiently  small,  one  can  drop  constants  B,a,(3, 
so  as  the  second  term  in  the  right  hand  side  of  (11.21). 

11.4  Stability  estimates  of  initial  conditions  of  parabolic  equations 
and  inequalities  in  infinite  domains 

The  result  of  this  subsection  was  obtained  by  the  PI  and  one  of  his  collaborators  Dr.  A.V. 
Tikhonravov  (Professor  of  Moscow  State  University)  [19].  Suppose  that  in  the  previous 
subsection  the  domain  12  C  M”  is  infinite  rather  than  finite,  <912  G  C 1 .  Consider  the  same 
problem  as  in  the  previous  subsection.  Let  the  function  u  G  H2'1  (Qt)  satisfies  the  following 
conditions 

Ut  =  Lu  +  / (x,  t),  a.e.  in  Qt,  (11.20) 

u|sT=0,  (11.21) 

u(x,  0)  =  g(x)  G  H1  (12) .  (11.22) 

Conditions  imposed  on  the  coefficients  of  the  elliptic  operator  L  are  the  same  as  ones  in 
the  previous  subsection.  Along  with  the  equation  (7.28)  consider  a  more  general  parabolic 
inequality 

| ut  -  L0u\  <  M  [|Vu]  +  \u\  +  |/|] ,  a.e.  in  QT ,  (11.23) 

where  Lq  is  the  principal  part  of  the  elliptic  operator  L  and  M  is  a  positive  constant.  Let 

P  G  C2,  P  C  12  be  an  arbitrary  surface  and  Pt  =  P  x  (0,  T).  Consider  the  following 
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Inverse  Problem.  Let  $  C  H  be  a  finite  subdomain.  For  either  of  problems  (11.20)- 
(11.22)  or  (11.21)-(11.23),  estimate  the  unknown  initial  condition  g(x )  in  the  domain  <f>  via 
the  lateral  Cauchy  data  h3(x,t)  and  h4(x,t),  where 

Qu 

U  I PT=  h3(x,t),—  I pT=  h4(x,t).  (11.24) 

The  main  difficulty  of  the  case  of  the  infinite  domain  12  compared  with  the  previous  case 
of  the  finite  domain  is  that  the  previously  used  idea  of  combining  lateral  and  backwards 
Carleman  estimates  does  not  work  in  this  case.  Indeed,  lateral  Carleman  estimates  work 
only  in  finite  subdomain.  Hence,  the  lateral  Carleman  estimate  would  enable  one  only  to 
estimate  the  norm  ||u(x,  £o)||l2(<j>)  for  a,  t0  =  const.  G  (0,  T) .  Then,  however,  one  cannot  use 
the  backwards  Carleman  estimate  (unlike  the  previous  subsection)  to  estimate  the  function 
g(x).  The  reason  of  the  latter  is  that  any  backwards  estimate  requires  the  knowledge  of  the 
function  u(x,to)  for  all  x  G  12.  This  is  why  analogues  of  Theorem  11.5  are  unknown. 

To  overcome  this  difficulty,  a  new  lateral  Carleman  estimate  for  the  parabolic  operator 
dt  —  Lq  was  derived.  The  level  surface  of  the  corresponding  Carleman  Weight  Function 
(CWF)  is  contained  in  a  thin  strip  t  €  {|f  —  <5|  <  ,  where  5  >  0  is  sufficiently  small 

and  the  number  ujq  G  (0, 1)  is  fixed.  The  main  new  feature  of  this  estimate  is  that,  unlike 
previously  known  Carleman  estimates,  this  one  does  not  break  down  when  the  width 
of  this  strip  approaches  zero  as  S  =  <5(||P||)  — >  0+  for  ||P||  — >  0. 

Before  formulating  the  stability  estimate,  we  formulate  a  geometric  condition.  Let  <f>  C  12 
be  a  convex  bounded  subdomain.  We  shall  say  that  <f>  has  the  P-property,  if  the  following 
two  conditions  are  fulfilled:  (1)  For  any  point  x  G  $  there  exists  a  point  x  (. x )  G  P  such  that 
the  straight  line  connecting  points  x  and  x  does  not  lie  in  the  hyperplane,  which  is  tangent  to 
the  hypersurface  P  at  the  point  x  and  (2)  dist  [<F,  (<912\P)]  >  0,  where  dist  [<f>,  (<912\P)]  :  = 
ds  (<f>)  is  the  Hausdorff  distance.  An  example  of  the  P-property  is  the  case  when  either 
P  C  (9<f>  or  P  C  <912  and  ds  ($)  >  0.  Another  example  is  when  the  hypersurface  P  is  a  part 
of  a  hyperplane,  P  C  {x\  =  0},  <f>  C  {x\  >  0}  D  H  or  <f>  C  [x\  <  0}  fl  O  and  ds  (<f>)  >  0.  The 
following  analogue  of  Theorem  11.4  was  proven 

Theorem  11.5  [19].  Suppose  that  above  conditions  imposed  on  coefficients  of  the  operator 
L(x,t,D),  the  domain  12  and  the  surface  P  are  fulfilled.  Let  the  function  u  G  H21  (Qt) 
satisfies  either  conditions  (11.20)-(11.22),  (11.24)  or  conditions  (11.21)-(11.24).  Let  <F  C  12 
be  a  convex  bounded  subdomain  of  the  domain  12  which  possesses  the  P-property.  Let  the 
function  h3  G  H 1  (Pr)  ■  Consider  the  vector  valued  function  F  =  (h3,  h4,  f )  and  denote 


11*11=  11*3 


1 2 

I  Hi  (Ft) 


+  II  Ml2(Pt)  + 


2 

^2  (Qt) 


1/2 


Suppose  that  ||F||  <  P,  where  B  is  a  positive  constant.  Choose  an  arbitrary  number  a  G 
(0,  2) .  Then  there  exist  constants  C  >  0  and  f3  G  (0, 1)  such  that  the  following  stability 
estimate  holds 


l2($)  ^ 


<k 


a 


llivslll 


1,2  (H) 


+ 


l2(  n\&) 


-i 

+  C 


(11.25) 
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Here  the  constant  C  =  C\  (L,  M,  <h,  P )  depends  on  the  operator  L,  the  constant  M  in  (6.31 ), 
the  domain  $  and  the  surface  P.  The  constant  (3  depends  on  the  same  parameters  as  ones 
listed  for  C ,  as  well  as  on  the  number  a. 

Hence,  if,  in  particular,  it  is  known  a  priori  that  the  function  g(x)  has  a  finite  support 
with  g{x)  =  0  for  x  G  H\<h,  then  the  term  ||fl,||i2(Q\$)  =  0  in  (11.25). 
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